Global optimality bounds for the placement of control valves in water supply networks

This manuscript investigates the problem of optimal placement of control valves in water supply networks, where the objective is to minimize average zone pressure. The problem formulation results in a nonconvex mixed integer nonlinear program (MINLP). Due to its complex mathematical structure, previous literature has solved this nonconvex MINLP using heuristics or local optimization methods, which do not provide guarantees on the global optimality of the computed valve configurations. In our approach, we implement a branch and bound method to obtain certified bounds on the optimality gap of the solutions. The algorithm relies on the solution of mixed integer linear programs, whose formulations include linear relaxations of the nonconvex hydraulic constraints. We investigate the implementation and performance of different linear relaxation schemes. In addition, a tailored domain reduction procedure is implemented to tighten the relaxations. The developed methods are evaluated using two benchmark water supply networks and an operational water supply network from the UK. The proposed approaches are shown to outperform state-of-the-art global optimization solvers for the considered benchmark water supply networks. The branch and bound algorithm converges to good quality feasible solutions in most instances, with bounds on the optimality gap that are comparable to the level of parameter uncertainty usually experienced in water supply network models.


Introduction
The efficient management of hydraulic pressure in pipes results in reduction of leakage (Lambert 2000;Wright et al. 2015) and risk of pipe failure (Lambert and Thornton 2011), and it is therefore one of the main operational challenges in water supply networks (WSNs). Here we consider pressure management using pressure control valves, which regulate pressure at their outlet. We investigate the problem of simultaneously optimizing the placement and operational settings of control valves in WSNs, where the objective is to minimize average zone pressure (AZP). AZP is used as a surrogate measure for leakage. The problem formulation includes flows across network links and hydraulic heads at nodes as continuous decision variables. In addition, binary variables are introduced to model the placement of valves. Mass and energy conservation laws are enforced as optimization constraints, resulting in a nonconvex mixed integer nonlinear program (MINLP). The solution of process network optimization problems frequently relies on the solution of MINLPs. Some examples include synthesis of heat exchanger networks (Zamora and Grossmann 1998), multi-period blending (Kolodziej et al. 2013), optimal design and operation of gas networks (Pfetsch et al. 2015;Humpola and Fügenschuh 2015), and water supply networks (DAmbrosio et al. 2015). In the framework of WSNs, MINLP formulations are ubiquitous and employed in a variety of applications, ranging from optimal network design (Bragalli et al. 2012;Sherali et al. 1999) to pump scheduling (Menke et al. 2015;Gleixner et al. 2012).
Both heuristic and mathematical optimization methods were applied in previous work to solve the problem of optimal valve placement in water networks. Heuristic approaches based on genetic algorithms (GAs) have been widely used for solving the considered problem-see Reis et al. (1997), Araujo et al. (2006), Nicolini and Zovatto (2009), Liberatore and Sechi (2009), Ali (2015), De Paola et al. (2017). However, they present some limitations. Firstly, they can not guarantee optimality of the computed solutions, not even local optimality. Moreover, the number of objective function evaluations and hydraulic simulations required by these approaches grows rapidly with the size of the network, precluding the application of GAs when large operational water network are considered. Previous work has also investigated the application of mathematical optimization methods for the solution of the problem of optimal valve placement in WSNssee Hindi and Hamam (1991), Eck and Mevissen (2012), Dai and Li (2014), Pecci et al. (2017a, b). Since the considered problem is nonconvex, approaches implemented in previous literature do not provide theoretical guarantees on the global optimality of the computed valve configuration. This paper investigates mathematical optimization methods to generate a certified bound on the optimality gap of the computed solutions for the problem of optimal valve placement in WSNs, guaranteeing -sub-optimality. We formulate a branch and bound method, which is a common approach in global optimization. To the best of the authors knowledge, global optimization methods have not been previously applied to the problem of optimal valve placement in WSNs. Previous literature has investigated global optimization techniques for optimal design of WSNs (Sherali et al. 1999;Raghunathan 2013). However, when pipe diameters are fixed, hydraulic heads and flow rates are uniquely determined and can be found by solving a strictly convex optimization problem (Raghunathan 2013). As a result, in the case of optimal WSN design, it is sufficient to focus the branch and bound on integer decision variables. On the contrary, when optimal operation of water supply networks is considered, spatial branch and bound is needed (Gleixner et al. 2012). In the present manuscript, we consider the problem of optimal valve placement in WSNs, where locations and operational settings of the control valves need to be simultaneously optimized. Therefore, branching is required on both continuous and integer variables.
The implemented branch and bound algorithm relies on a sequence of lower and upper bounds to the optimal value of the nonconvex MINLP in study-for a general review see Tawarmalani and Sahinidis (2002). Since all convex constraints within the problem formulation for optimal valve placement are linear, it is particularly convenient to generate lower bounds using linear relaxations of the nonconvex constraints (Tawarmalani and Sahinidis 2002, Chapter 4). The nonconvexity of MINLPs arising in the framework of water networks is due to the absolute power functions representing friction energy losses within the system's conservation laws-see Eq. (2a). Analogous nonconvex expressions have previously been studied in other engineering frameworks, where linear relaxations were formulated-see Humpola and Fügenschuh (2015), Gleixner et al. (2012), Liberti and Pantelides (2003), Tawarmalani and Sahinidis (2002), Udell and Boyd (2015), Vigerske (2012). We define linear relaxations of the nonconvex equality constraints considered here by extending the formulation proposed in Liberti and Pantelides (2003) for monomials of odd degree. Such linear relaxations define an outer approximation of the convex envelopes of the nonconvex equality constraints. We investigate the use of different number of linearizations for the outer approximation. The strength of the linear relaxations depends on the diameter of the decision variables' domain (Puranik and Sahinidis 2017). Therefore, we implement a domain reduction procedure, based on the solution of a series of linear programs (LPs). The proposed approach takes advantage of the underlying network structure to reduce the number of linear programming solves. Benefits and limitations of the developed methods are evaluated using two benchmark water networks, and a large-scale operational network from the UK. Moreover, the numerical results show that the proposed approach outperforms state-of-the-art global optimization solvers for the considered benchmark water networks. The branch and bound framework has enabled the convergence to -sub-optimal solutions for the problem of optimal valve placement, with bounds on the optimality gap comparable to the order of parameter and data uncertainties inherent in operational network models.

Problem formulation
A water supply network with n 0 water sources, n n demand nodes and n p pipes, is modelled as a directed graph with n n + n 0 nodes and n p edges. The operation of a network is considered under n l different demand conditions during the diurnal cycle.

3
The nodal demands are denoted by d t ∈ ℝ n n , while known hydraulic heads at water sources are indicated by h t 0 ∈ ℝ n 0 , for each t = 1, … , n l . Furthermore, the vector of node elevations is represented by ∈ ℝ n n . Given t ∈ {1, … , n l } , we consider hydraulic heads h t ∈ ℝ n n and flow rates q t ∈ ℝ n p as continuous decision variables. Moreover, vector t ∈ ℝ n p is included to model the unknown head loss introduced by the action of pressure control valves. We introduce auxiliary variables t ∈ ℝ n p to isolate the nonconvex terms representing the friction head losses occurring within the pipes of a network. These can be expressed by either the Hazen-Williams (HW) or Darcy-Weisbach (DW) formulae. Since both friction head loss formulae involve non-smooth nonconvex terms, it is convenient to use smooth quadratic approximations, computed over a range of flow (Pecci et al. 2017c). When a suitable quadratic approximation has been determined, it can be written as j (q t j ) = q t j (a j |q t j | + b j ) , with a j ≥ 0 and b j ≥ 0 , for all j = 1, … , n p and t = 1, … , n l -see also Eq. (28) in Appendix 1.
In this work, we investigate the optimal placement and operation of pressure control valves to minimize average zone pressure (AZP), defined as: L j 2 and I(i) is the set of indices for links incident at node i, L j is the length of pipe j and W = ∑ n n i=1 w i is a normalisation factor. The optimization of valve placement and operation to minimize AZP is primarily subject to conservation energy (2b) and mass (2c) laws: where matrices A 12 ∈ ℝ n p ×n n and A 10 ∈ ℝ n p ×n 0 are edge-node incidence matrices for unknown and known head nodes, respectively. Moreover, we define the vector of friction head losses as We consider vectors of binary variables z + , z − ∈ {0, 1} n p with the following properties: prevents the placement of two valves on the same link.
The placement of control valves in a WSN is modelled by big-M constraints, ensuring that additional head losses have the same direction as flows and friction head losses across the valves. The following constant vectors and matrices are used to formulate the big-M constraints. Given a vector of maximum allowed velocities across network pipes v max ∈ ℝ n p , let q t L = −Sv max , and q t . Minimum and maximum allowed hydraulic heads at nodes are specified by h t min , h t max ∈ ℝ n n , respectively. Define t L ∈ ℝ n p ad t U ∈ ℝ n p as follows: . We formulate the following constraints: Let n v be the number of valves to be installed. The desired properties of the binary variables are enforced by the linear constraints: Finally, physical and operational bounds on all continuous variables are enforced by the following constraints: (3)

3
The valve placement optimization problem can be written in a compact form as: is a mixed integer nonlinear program (MINLP). Furthermore, the equality constraints (2a) are nonconvex, hence the problem formulation results in a nonconvex MINLP.
We implement a branch and bound method that relies on the generation of a sequence of lower and upper bounds to the optimal value. The present work takes a similar approach to Floudas (2013, 2014) and compute lower bounds solving Mixed Integer Linear Programming (MILP) relaxations of MINLP(Q) . As discussed in Smith and Pantelides (1999), MILP relaxations result in tighter lower bounds then relaxed linear programs (LPs) and the algorithm is expected to converge in less iterations. However, they also require an higher computational effort for each branch and bound iteration. Nonetheless, as shown in the numerical results reported in Sect. 4, the implementation of state-of-the-art MILP solvers (e.g. Gurobi Optimization 2017) have enabled the application of the considered method to large problem instances.

Lower bounding MILP
for all t ∈ {1, … , n l } . We consider the following restriction of MINLP(Q): where F(Q � ) ⊂ ℝ n l n p × ℝ n l n n × ℝ n l n p × ℝ n l n p × {0, 1} n p × {0, 1} n p is the set defined by (2b)-(6d) using the bounds in Q ′ . Let y * (Q � ) be the optimal value of MINLP(Q � ) . In the following, we formulate a MILP relaxation of MINLP(Q � ) , which yields a lower bound L(Q � ) ≤ y * (Q � ).
A detailed derivation of linear relaxations of (2a) is presented in Appendix 1, where the formulation for monomials of odd degree proposed by Liberti and Pantelides (2003) is extended to functions ( j (⋅)) j∈{1,…,n p } . Such linear relaxations can be written as for suitable matrices R t and E t , and vectors r t , which depend on (q t L ) � , (q t U ) � , and a parameter N c ≥ 0 used to control the number of linearizations for the outer approximation of the convex envelopes of constraints (2a)-see Appendix 1. In Section 4, we investigate the effect of different number of linearizations on the convergence properties of the branch and bound algorithm. Using the above linear relaxations, it is possible to define the following MILP relaxation of MINLP(Q � ): Let L(Q � ) be the optimal value of MILP(Q � ) . If MILP(Q � ) is infeasible, then MINLP(Q � ) is infeasible and we have L(Q � ) = y * (Q � ) = +∞ . Otherwise, the solution of MILP(Q � ) yields a lower bound L(Q � ) ≤ y * (Q � ) . Moreover, if the current best lower bound LB is available, set L(Q � ) ← max(L(Q � ), LB).
The generation of a valid upper bound via the (local) solution of a nonlinear program is often computationally expensive, particularly when large problem instances are considered. However, NLP(Q � ) presents a high level of sparsity, that is retained by constraints (2b) and (2c) from the sparse structure of water supply networks. As a consequence, sparse NLP solvers offer scalable solution approaches for NLP(Q � )-see the interior point method introduced in Waechter and Biegler (2006).

Domain reduction
The strength of the MILP relaxations is expected to have a significant impact on the convergence properties of branch and bound schemes (Belotti et al. 2009). This is confirmed by the numerical results reported in Sect. 4. As discussed in Appendix 1, in the case considered here, smaller ranges of flows lead to tighter relaxations. Therefore, we investigate pre-processing strategies to reduce the domains of the flow variables, focusing on Optimization Based Bound Tightening (OBBT), which relies on the solution of a series of optimization problems to tighten upper and lower bounds on selected variables-for a recent review on variable bound tightening approaches see Puranik and Sahinidis (2017).
In order to highlight the separable structure of the feasible set of MINLP(Q) with respect to the time indices, we introduce the time indexed sets {G t (Q)} t=1,…,n l such that G t (Q) ∈ ℝ n p × ℝ n n × ℝ n p × ℝ n p × {0, 1} n p × {0, 1} n p is defined by where z + = u t 1 and z − = u t 2 , for all t ∈ {1, … , n l } , and q = (q t ) t=1,…,n l , h = (h t ) t=1,…,n l , = ( t ) t=1,…,n l , = ( t ) t=1,…,n l . By definition of {G t (Q)} t=1,…,n l , the feasible set of MINLP(Q) is equivalent to the feasible set of the following problem: Global optimality bounds for the placement of control valves… Let ∈ {−1, 1} , t ∈ {1, … , n l } , and j ∈ {1, … , n p } . Consider the mixed integer linear program: where R l , E l and r l represent the linear relaxations of (2a) computed using the flow bounds in Q as described in Appendix 1. The feasible set of Problem (10) is a relaxation of the set of feasible solutions of Problem (9). Tightening upper and lower bounds on the flow variables would require the solution of 2n l n P mixed integer linear programs analogous to (10). Mixed integer linear programs can be solved by state-of-the-art MILP solvers (Gurobi Optimization 2017). However, the dimension of Problem (10) and the number of problems to be solved can pose significant computational challenges when considering large-scale water networks. Therefore, the remainder of this section investigates tailored strategies to reduce the computational effort required by the considered domain reduction procedure. Firstly, Problem (10) is converted into a linear program (LP) by replacing G t (Q) with a valid polyhedral relaxation Ĝt (Q) . Furthermore, the resulting LP can be decoupled with respect to the time indices l ∈ {1, … , n l } , omitting the consistency constraints in Problem (10). For each t ∈ {1, … , n l } and j ∈ {1, … , n p } , we consider the following relaxation of Problem (10): The OBBT method solves LP If Q ′ is the resulting set of bounds, then MINLP(Q � ) and MINLP(Q) have the same optimal solution. Further domain reductions are achieved by iteratively applying the described process. In order to reduce the number of flow variables whose bounds are tightened through the solution of linear programs, we exploit the underlying network (graph) structure of LP ,t,j Q . In what follows, few graph-theoretical definitions for WSNs are reviewed. The degree of an unknown head node in a WSN graph is defined as the number of links connected to the node. We define a tree in a WSN graph as an acyclic connected subgraph such that at least one of its unknown head nodes has degree one, and only one of its nodes is connected to either a looped part of the network or to a fixed head node (Deuerlein 2008). Such a unique node of a given tree is called its root node. The forest of a water network is 1 3 defined as the disjoint union of all trees in the network (Deuerlein 2008). The part of the network which is not contained in the forest but includes the roots of all the trees is called core. If a link belongs to the forest of a network graph, its flow is uniquely determined by the demand assigned to the adjacent forest nodes (Elhay et al. 2014). As a consequence, flows across forest links do not depend from the optimization process and their values are fixed. The forest-core decomposition of a WSN graph can be implemented using the algorithm presented in Elhay et al. (2014).
In the bound tightening approach proposed here, upper and lower bounds on flow variables corresponding to forest links are set a priori. Then, let core links j 1 and j 2 be connected in series through a node i. For each t ∈ {1, … , n l } , flow variables in LP ,t,j Q are subject to the following mass conservation laws: Hence: Given a sequence of core links connected in series, it is possible to select one link as representative and derive relations analogous to (12) for all other links in the sequence, proceeding by substitution. Let ⊂ {1, … , n p } be a subset of indices corresponding to network's links belonging to the core, where a unique representative for each sequence of links have been selected. An illustrative example of the definition of for a simple network model is included in Appendix 2.
The proposed strategy is to perform OBBT only for those flow variables whose indices are in and then propagate the bounds to the remaining core links using relations analogous to (12). Algorithm 1 is implemented for iterative domain reduction on the flow variables in MINLP(Q) . Such iterative method can fail to converge to a fixed point in finite time (Puranik and Sahinidis 2017). Therefore, Algorithm 1 is terminated when the progress in the domain reduction is not significant or when the maximum number of iterations is reached. Given a possible choice of lower and The analytical study of the number of iterations required to generate tight variable bounds represents an open research problem (Puranik and Sahinidis 2017). Therefore, the maximum number of iterations in Algorithm 1 is set to 10 based on empirical experience. The number of linear programming solves at each iteration of Algorithm 1 is 2n l | | , where | | is expected to be significantly smaller than n p when operational water networks with a relatively low number of loops are considered. Furthermore, the for-loops within Algorithm 1 can be executed in parallel, as each iteration does not depend from the others.

Branch and bound algorithm
Let Q init be the set of initial lower and upper bounds on flow variables. It can either be Q init = Q or Q init = Q tight , where Q tight represents tightened lower and upper bounds on flow variables resulting from Algorithm 1. We denote by y * = y * (Q) the optimal value of MINLP(Q) . Then, by definition, y * = y * (Q init ) . The branch and bound method iteratively generates a hierarchy of optimization problems, starting from MINLP(Q init ) , represented by a binary tree named branch and bound tree. In particular, the algorithm creates a partition of F(Q init ) given by {F(Q � ) | Q � ∈ } such that As standard in branch and bound methods, the proposed algorithm starts by computing lower and upper bounds to y * -this is iteration 0. Then, the root problem MINLP(Q init ) is branched into two optimization problems-see Sect. 3.5. After branching, the algorithm computes lower and upper bounds on the two descendant problems as outlined in Sects. 3.1 and 3.2, respectively. The best upper and lower bounds are updated as necessary. Finally, a new problem from the branch and bound tree is selected for branching and the iteration is repeated. The method stops when a termination criterion is met and its implementation is described in Algorithm 2. The output of the algorithm is either an -sub-optimal solution to MINLP(Q) or a feasible solution with a certified bound on the optimality gap greater than .
The computed optimality gaps should be considered within the range of uncertainties that are inherent in modelling of operational water networks. As discussed in Wright et al. (2015), pressure control in operational water networks is subject to multiple sources of data and modelling errors. These include the stochastic nature of customer demand, uncertainty in the hydraulic model parameters, network connectivity, measurements accuracy, and factors affecting the physical operation of control and isolation valves.
Remark 1 If Q � ∈ is such that L(Q � ) > UB , then the global optimum cannot belong to F(Q � ) and MINLP(Q � ) can be pruned from the branch and bound tree. However, the selection strategy implemented in Algorithm 2 implies that the branch and bound iterations will terminate before Q ′ is selected for branching. Therefore, the method does not explicitly implement a pruning procedure.

Branching strategy
Let MINLP(Q b ) be a problem in the branch and bound tree such that . A good branching strategy should aim to improve the lower bound (i.e. L(Q left ) ≥ L(Q b ) and L(Q right ) ≥ L(Q b ) ), and maintain a balanced branch and bound tree (Belotti 2013, Section 5.3.1). A branch and bound tree is said balanced if all problems in the tree are similarly difficult to solve. In this work, the following branching rule is implemented.
,n l be as follows: Then, define of F(Q left ) and F(Q right ) using the bounds on flow variables contained in Q left and Q right .
The above branching strategy improves the lower bound, by definition of the linear relaxations given in Appendix 1. In fact, the formulation of MILP(Q left ) and MILP(Q right ) requires the inclusion of linear inequalities corresponding to the polyhedral relaxations computed with the new bounds on variable q l k . Such refined linear relaxations are exact at q l k resulting in tighter MILP relaxations. If q l k is too close to either (q l L ) b k or (q l U ) b k , the implemented branching strategy can result in an unbalanced branch and bound tree. However, the relaxation error Fig. 11. As a result, in most cases, q l k is expected to be closer to the middle point of the interval

Results and discussion
In this section, the developed methods are evaluated using two benchmark network models, and a large operational water supply network from the UK. The corresponding problem sizes are reported in Table 1. All experiments presented below were conducted within MATLAB 2016b-64 bit for Windows 7, installed on a 2.50 GHz Intel Xeon(R) CPU E5-2640 0 with 12 Cores and 12 GB of RAM. The MILPs and LPs involved in Algorithms 1 and 2 are solved using GUROBI (v7.5) (Gurobi Optimization 2017), which is accessed via its MATLAB interface. In addition, all iterations in Algorithm 1 were executed in series, and GUROBI was forced to operate on a single-thread; all other parameters in GUROBI are set to their default values. Local solutions to the upper bounding NLPs considered in Sect. 3.2 are computed using IPOPT (v.3.12.6) (Waechter and Biegler 2006), accessed in MATLAB via an interface of the OPTI Toolbox (Currie and Wilson 2012). In the implementation of IPOPT, sparse gradients and Jacobians are directly supplied to the solver, in order to take advantage of the very sparse structure of our problem. The global optimization solvers BARON (v18.8.23) (Klnç and Sahinidis 2018) and SCIP (v3.2.1) (Gamrath et al. 2016) were implemented for the direct solution of the considered nonconvex MINLPs. The solver SCIP is accessed in MATLAB via an interface provided by OPTI Toolbox and relies on SoPLEX (v221) as linear solver and IPOPT (v3.12.6) as nonlinear solver. In our implementation, we accessed BARON via its MATLAB interface, and used IBM ILOG CPLEX (v12.8) as linear solver, while BARON was allowed to select the nonlinear solver according to its dynamic strategy (Sahinidis 2018). Finally, we set = 10 −6 , MaxIter = +∞ , and define ∶= 100 ⋅ UB−LB LB . Firstly, we consider PescaraNet, a reduced model of the water supply network for the city of Pescara (Bragalli et al. 2012). The network's layout is presented in Fig. 1a. This network model has 68 demand nodes, 99 pipes and 3 water sources. Moreover, friction head losses are modelled using the HW formula. In this case, the formulation of MINLP(Q) includes a single demand condition and results in a relatively small size nonconvex MINLP whose characteristics are reported in Table 1. This water network has been already considered in Eck and Mevissen (2013) where a local optimization method was applied for computing optimal valve locations in PescaraNet. However, Eck and Mevissen (2013) did not include details on upper and lower bounds of flow variables, which precludes a direct comparison with the solutions obtained here. In the present implementation, maximum allowed velocity in each pipe is set to 2 (m/s) and minimum required pressure at each node to 19 (m). For every link j ∈ {1 … , n p } , a quadratic head loss approximation j (⋅) is defined as in Pecci et al. (2017c).
The second case study is named Net25. This benchmark network model has been used to evaluate solution approaches for optimal valve placement problems in previous literature using heuristics (Reis et al. 1997;Araujo et al. 2006;Liberatore and Sechi 2009;Nicolini and Zovatto 2009;Ali 2015;De Paola et al. 2017) and mathematical optimization methods (Eck and Mevissen 2012;Dai and Li 2014;Pecci et al. 2017a, b). The results presented in this section allow the quantification of the level of sub-optimality of valve configurations for Net25 previously computed using heuristics and local optimization methods. The network has 22 nodes, 37 pipes and 3 reservoirs-see Fig. 1b. Details on pipes' characteristics, nodal demands and reservoirs' levels are presented in Jowitt and Xu (1990) and Dai and Li (2014). The hydraulic model of Net25 uses the HW formula to model friction head losses. We observe that Net25 has a smaller dimension with respect to PescaraNet; nonetheless, it results in a larger nonconvex MINLP as the problem formulation considers 24 demand conditions, one for each hour of the day-see Table 1. Moreover, we set the maximum allowed velocity in each pipe to 1 (m/s) and the minimum pressure at each node to 30 (m). Analogously to what done for PescaraNet, for each link j, a quadratic approximation of the friction losses j (⋅) is computed as proposed in Pecci et al. (2017c).
The developed methods are tested for optimizing the placement of 1 to 5 control valves in both PescaraNet and Net25, where the objective to be minimized is AZP.
In these experiments, we set a time limit of 7200 (s). Given > 0 , let ( ) be the percentage of problems with remaining smaller or equal than , defined as

3
For the purpose of this numerical study, we consider a uniform distribution of values of between 0 and 100, spaced by 10 −4 . Let N c ≥ 0 be the parameter controlling the number of linearizations for the outer approximation of the convex envelopes of constraints (2a) -see Appendix 1. In the following, the symbol BB N c indicates the implementation of Algorithm 2 where the polyhedral relaxations (7) are computed as detailed in Appendix 1. When Algorithm 1 is implemented as pre-processing routine to reduce the domain of the flow variables, the combination of Algorithms 1 and 2 is denoted by dBB N c . Moreover, in order to ensure a total maximum computational time of 7200 (s), the value of T max in Algorithm 2 is adjusted to account for the time spent in Algorithm 1. Initially, we evaluate the performance of BB N c with N c ∈ {0, 1, 3, 5} . The complete numerical results are reported in Tables 2, 3, 4, 5, 6, 7, 8, and 9 of Appendix 3. When no domain reduction is applied, the solutions computed after 7200 (s) by BB N c are within 10% of optimality, for all N c ∈ {0, 1, 3, 5} . However, in the case of Net25, the obtained valve placements for n v ∈ {3, 4, 5} differ from the best-known solutions for the case study, with respect to those reported in Eck and Mevissen (2012), Dai and Li (2014), Pecci et al. (2017a), Pecci et al. (2017b). As shown in Fig. 2a, the inclusion of more linearizations for the outer approximation of the convex envelope of (2a) improves the convergence properties of the branch and bound algorithm with respect to BB 0 . In addition, we observe that algorithms BB 1 , BB 3 , and BB 5 had similar performances. However, only BB 1 terminated in all instances with a Gap smaller than 7%.
We investigate the implementation of algorithms dBB N c , with N c ∈ {0, 1, 3, 5} , for solving the considered optimal valve placement problems in PescaraNet and Net25. Firstly, Algorithm 1 was applied to tighten the flow variable bounds, for each  Tables 10, 11 , 12, 13, 14, 15, 16, 17. As shown in Fig. 2b, in most instances, the solutions computed by dBB N c are within 5% of optimality, for all N c ∈ {0, 1, 3, 5} . Moreover, all the computed feasible solutions for Net25 equal the best-known valve locations obtained for the considered case study with respect to those reported in Eck and Mevissen (2012), Dai and Li (2014), Pecci et al. (2017a, b). According to these numerical results, the implementations of dBB N c with N c > 0 have better computational performance than dBB 0 . In order to further investigate the effect of the domain reduction procedure on the convergence properties of Algorithm 2, we compare the implementation of the branch and bound algorithm with and without tightening the bounds on flow variables. As shown in Figs. 3, the domain reduction procedure significantly improved the performance of the branch and bound method for the considered problems of optimal valve placement in PescaraNet and Net25, irrespectively of the number of linearizations used. Furthermore, as shown in Fig. 4, the implementation of dBB 0 achieves better performance than the direct applications of branch and bound algorithm with N c ∈ {1, 3, 5} . These results suggest that, for these benchmark water networks, increasing the number of linearizations improves the branch and bound ability to reduce the optimality gap, but not as much as the implementation of the developed domain reduction procedure. Furthermore, all the numerical experiments show rapid convergence of the algorithms to the final feasible solution and the associated optimality gap. As examples, Fig. 5 reports the progress of dBB 0 and dBB 1 when solving the problems of optimal placement of 3 valves on PescaraNet and Net25, respectively. In particular, it shows that good quality (if not optimal) solutions are computed early in the algorithm, while the remaining iterations are mostly used to improve the lower bounds.
In conclusion, we compare the performance of dBB 5 with state-of-the-art global optimization solvers SCIP and BARON for the solution of the considered nonconvex MINLPs. A time limit of 7200 (s) was set for all solvers. All other options in BARON and SCIP are set to their default values. As reported in  Tables 20 and 21 show that SCIP has failed to find a feasible solution only for the problem of optimal placement of 5 valves in Net25. However, in most instances, SCIP resulted in a larger than 5% . Furthermore, the majority of the feasible solutions computed by SCIP differ from the bestknown solutions for Net25, which are reported in Pecci et al. (2017a) and were computed by the convex MINLP solver BONMIN (Bonami et al. 2008). Results reported in Tables 16, 17 , 18, 19, 20, 21 show that the lower bounds computed by dBB 5 are tighter than those obtained by BARON and SCIP, in all instances. As illustrated in Fig. 6, after two hours of computations, the developed branch and bound method consistently results in smaller optimality gaps than the two stateof-the-art solvers.

Case study 3: BWFLnet
Finally, the proposed global optimization method for optimal valve placement has been applied to BWFLnet, the network model of the Smart Water Network Demonstrator (Field Lab) operated by Bristol Water, InfraSense Labs at Imperial College London and Cla-Val (Wright et al. 2014). This water supply network consists of 2310 nodes, 2369 pipes and 2 inlets (with fixed known hydraulic heads); its graph is presented in Fig. 7. The HW formula is used to model friction losses within BWFLnet. The size of BWFLnet is two orders of magnitude larger than the other considered benchmark WSNs. In BWFLnet, the network's operator has already installed three automatic pressure control valves and two boundary control valves, which are optimally controlled in order to minimize AZP (Wright et al. 2015). In the problem formulation for optimal valve placement in BWFLnet, we model the three pressure control valves as open smooth pipes. Moreover, details on the operation of the two boundary valves have been provided by the valves' manufacturer. The daily demand profiles in BWFLnet include 96 different demand conditions, one every 15 minutes. However, the problem of optimal valve placement and operation in BWFLnet under 96 different demand conditions would result in a large MINLP with 227, 424 nonconvex constraints. In order to limit the problem size, MINLP(Q) is formulated on BWFLnet using only three most relevant demand conditions, which were selected to capture network's typical daily operational conditions at low night time demand, morning peak demand, and afternoon dip-see Fig. 8a. The AZP values computed using a problem formulation restricted to these 3 demand conditions are expected to be close to those obtained for a full set of 96 different demand conditions-this is confirmed to be the case for BWFLnet, as shown in Tables 24 and 26.
Note that the highest impact in terms of AZP reduction will be achieved by controlling pressure through pipes carrying large quantity of water, hence links with small diameters are not likely to be good candidates for control valve placement. As a consequence, links in BWFLnet whose diameter is smaller than 100 (mm) are discarded from the set of candidate valve locations. The minimum allowed pressure at all demand nodes is set to 18 (m), while this value is relaxed to zero for nodes with no demand. A quadratic approximation of friction losses is computed as discussed in Pecci et al. (2017c). The size of the resulting MINLP is shown in Table 1 b Maximum velocity achieved in simulation across each link expected to be significantly smaller than 3 (m/s). As a consequence, when upper and lower bounds on the flow variables are based on a maximum allowed velocity of 3 (m/s) for all pipes, the polyhedral relaxations (7) are not expected to be tight, leading to loose MILP relaxations-see also the discussion in Appendix 1. In order to tighten the polyhedral relaxations (7) and avoid unnecessarily large flow bounds, we define a tailored maximum allowed velocity for each link. It is known that the placement of control valves can result in velocity changes across network links (Abraham et al. 2018). In order to limit the possibility of discarding optimal solutions from the feasible set, we implement the following heuristic. Given a link j ∈ {1, … , n p } , let sim j be the maximum velocity achieved in simulation under multiple demand conditions. If sim j ≤ 1 ( m s ) , then the optimized velocities across link j are allowed to exceed sim j by at least a factor of two. In comparison, for links where the simulated velocities achieved higher rates, the allowed increment is more limited. The implemented strategy is detailed in the following: Then, MINLP(Q) was formulated on BWFLnet using tailored maximum velocities to define upper and lower bounds on flow variables Q. We consider the optimal placement of 1 to 5 control valves in BWFLnet. In the numerical experiments reported in this section, we set a time limit T max = 86,400 (s) (1 day). We investigate the direct application of Algorithm 2 with N c ∈ {0, 3} , without implementing the domain reduction procedure. The results are summarised in Tables 22 and 23. In the case of N c = 0 , the relative optimality gap obtained for n v ∈ {1, 2, 3} are larger than 30% . No feasible solution was found after a day of computations for n v ∈ {4, 5} . Including more linearizations for the outer approximation of the convex envelopes of (2a) results in small improvements for the cases of n v ∈ {1, 2} , while no feasible solution was obtained when n v ∈ {3, 4, 5}-see Table 23. Then, the domain reduction procedure described in Algorithm 1 was applied as pre-processing reoutine to tighten the flow variable bounds within the formulation of MINLP(Q) on BWFLnet. In this case, the fraction of links involved in LP solves corresponds to less than the 10% of the whole set of network links. Such reduction is due to the structure of BWFLnet, a typical water network from the United Kingdom, where a considerable portion of links belongs to the forest. Moreover, hydraulic models of operational networks often present sequences of links connected in series, used to model the existence of different customer connections along network pipes. As a result, the proposed graph-theory based decomposition has significantly reduced the computational cost associated with Algorithm 1, which required the solution of 2532 LPs for each case of n v = 1, … , 5 . In the numerical experiments reported here, such LPs were solved sequentially. As a result, Algorithm 1 required roughly one hour of CPU time, in all instance-see Tables 24 and 25. However, as previously observed, the LPs solved at each iteration of Algorithm 1 do not depend on each other. Therefore, in a practical implementation, they can be solved in parallel, exploiting the existence of multiple computational cores.
Algorithm 1 was then applied to MINLP(Q tight ) , with N c ∈ {0, 3}-see Tables 24  and 25. When N c = 0 , the bounds on the absolute optimality gaps obtained for n v ∈ {1, 2, 3} are between 3 and 5 (m)-see also Fig. 9. Optimality bounds of such magnitude are comparable to the order of uncertainty affecting pressure control in BWFLnet (Wright et al. 2015). Hence, the quality of the computed feasible solutions is considered to be acceptable. Again, in the cases of n v ∈ {4, 5} , no feasible solution was found before reaching the time limit. Similar optimality gaps were obtained when N c = 3 , but the algorithm failed to compute a feasible solution for n v ∈ {3, 4, 5} . These results suggest that the inclusion of different number of linearizations for the outer approximation of the convex envelopes of (2a) did not improve the performance of the branch and bound algorithm for solving the problem of optimal valve placement in BWFLnet.
We observe that the feasible solutions computed with and without domain reduction are the same. However, as shown in Fig. 9, the application of Algorithm 1 has significantly improved the quality of the lower bounds computed by the branch and bound algorithm. Nonetheless, the optimality gaps reported in Table 24 are too wide to be reduced within the prescribed time limit, even if the domain reduction routine is applied. The numerical results show that, in all problem instances, most pipes in BWFLnet experience low velocities, which are significantly smaller than their expected maximum velocities, despite the implemented heuristic to tailor maximum allowed velocities. As a consequence, the polyhedral relaxations of Eq. (2a) corresponding to these pipes are not sufficiently tight-see also the discussion in Appendix 1.
Next, given optimal valve locations for n v ∈ {1, 2, 3} , computed with Algorithm 2, we have formulated a nonconvex nonlinear program to optimize valve operation under the complete set of 96 demand conditions. Such NLP is obtained from the formulation of MINLP(Q) by fixing the values of the binary variables corresponding to optimized valve locations, and it is equivalent to the problem of optimizing the actuators operation. Algorithm 1 is applied to tighten the bounds on flow variables, and Algorithm 2 is subsequently implemented to solve the nonconvex NLP with tightened variable bounds. Lower and upper bounds computed at the root of the branch and bound tree are reported in Table 26, for n v ∈ {1, 2, 3} . As expected, the upper bounds to the AZP values computed considering the full set of multiple demand conditions are close to those obtained for the restriction to three demand conditions. Moreover, in this case, the domain reduction procedure has resulted in considerably tighter lower bounds, without performing any branching operation. In contrast, the results reported in Table 24 indicate that that inclusion of valve locations as unknowns significantly reduces the ability of Algorithm 1 to compute tight estimates of flow variable domains. This situation has an interpretation from the hydraulic application perspective. In fact, appropriate changes in network topology induced by valves closures can result in increased flow velocities across some pipes (Abraham et al. 2018).
In conclusion, we observe that the average computational cost associated with the solution of a single MILP relaxation for BWFLnet is considerably higher than for Net25 and PescaraNet. In fact, solving the MILP relaxation at the root of the branch and bound tree for BWFLnet required 2-6 orders of magnitude more computational time than what experienced for Net25 and PescaraNet. This behaviour is predictable, as the problem of optimal valve placement and operation in BWFLnet results in a MINLP whose size is one to two orders of magnitude larger than the size of the problem formulation on Net25 and PescaraNet-see Table 1. It is known that computational effort required to solve a mixed-integer program grows combinatorially with the size of the problem. These numerical experiments were conducted on a single computational thread of a desktop machine, using a standard implementation of GUROBI for the solution of the MILP relaxations. The availability of additional computational capability and the use of a tailored MILP solver could speed-up the optimization process.
All three case studies show that good quality solutions are often computed at the root node of the branch and bound tree (i.e. iteration 0). This is in accordance with the work by Diamond et al. (2018), presenting a number of examples where good quality solutions to nonconvex optimization problems are recovered from the solution of suitable convex relaxations. The results suggest that Algorithm 2 can be early terminated to generate good quality solutions of MINLP(Q) , by opportunely setting time limit or maximum number of iterations. Moreover, observe that Algorithm 2 provides more information than local optimization methods for optimal valve placement in WSNs studied in previous work. In fact, the algorithm always generates a certified bound on the optimality gap of the computed solution. Local optimization methods can be implemented before starting Algorithm 2, in order to rapidly generate good quality feasible solutions.

Conclusions and future work
In this manuscript, we have investigated the application of branch and bound strategies to compute -sub-optimal solutions for the problem of optimal valve placement and operation in water supply networks. The implemented algorithm relies on the solution of MILP relaxations of the original nonconvex MINLP. Furthermore, a tailored domain reduction procedure was implemented to tighten the MILP relaxations. In contrast to previously published solution methods for optimal valve placement in water networks, the presented algorithm terminates with a certified bound on the optimality gap of the computed solution, thus providing additional information to support the design and operation of water supply networks.
The proposed branch and bound method has successfully generated good quality feasible solutions in two benchmark water networks and a large water supply network, after few iterations, with bounds on the optimality gap comparable to the order of uncertainty usually experienced in pressure control of water supply networks. Furthermore, the results suggest that the proposed domain reduction strategy is more effective in improving the convergence properties of the algorithm, than simply increasing the number of linearizations used to define the polyhedral relaxations. Moreover, the results reported in this manuscript show that, for the considered benchmark water networks, the proposed branch and bound algorithms outperform state-of-the-art global optimization solvers SCIP and BARON. These results highlight the challenges of applying off-the-shelf global solvers to the problem in study.
However, the lower bounds generated by the algorithm experience slow progress, as shown in Fig. 5. They can be improved by performing Algorithm 1 on Q left and Q right , so that the new bounds on the branching variable can be propagated to the remaining flow variables. Recall that, at each iteration of Algorithm 1, the required solutions of the 2n l | | linear programs can be computed in parallel. As a consequence, if enough computational cores are available, the outlined bound propagation strategy can be applied at each branch and bound iteration, without dramatically increasing the overall computational time. Furthermore, future work should investigate the inclusion of additional valid linear inequalities in the formulation of MILP(Q � ) . A possible strategy is to use the locally optimal solutions generated at each stage of the branch and bound method, following an approach similar to Humpola et al. (2016). In addition, when large operational water networks are considered, tailored solvers for the relaxed MILPs can be designed using suitable decomposition strategies, following ideas discussed in Vigerske (2012, Chapters 3 and 4).
Although the present work has focused on the minimization of AZP, other convex objective functions can be minimized within the same framework, with little modifications to the discussed formulation. In particular, if the objective function is nonlinear, the generation of lower bounds requires the solution convex MINLPs. Furthermore, optimal design (Bragalli et al. 2012), optimal valve control (Wright et al. 2015), and optimal pump scheduling (Menke et al. 2015) problems result in optimization problems involving nonconvex constraints like (2a). As a consequence, using suitable linear relaxations and NLP sub-problems, Algorithms 1 and 2 can be applied to guarantee the optimality of the solution found, or a bound on the optimality gap, for such nonconvex optimization problems.
When q ≤ q L , it is necessary to modify the definition of (⋅) as follows (see Fig. 10b): Analogously, if q ≥ q U set (see Fig. 10c) Now, assume that 0 ≤ q L < q U . In this case, the convex envelope (see Fig. 10d) is given by (36) (q) = (q)

Fig. 10 Convex envelopes of nonconvex constraints
Analogously, if q L < q U ≤ 0 , set (see Fig. 10e) In this manuscript, our aim is to apply a branch and bound method to the nonconvex problem MINLP(Q) . At each step of the algorithm, a solution of a convex MINLP relaxation of MINLP(Q) will yield a lower bound on the optimal objective function value. Observe that, with the exception of (2a), all the other constraints in MINLP(Q) are linear. Therefore, it is convenient to study polyhedral relaxations of the nonconvex set in (26), which will result in a Mixed Integer Linear Programming (MILP) relaxation of MINLP(Q)-for a general review, see Tawarmalani and Sahinidis (2002, Chapter 4). Let N c ≥ 0 be a parameter used to control the number of linearizations for the outer approximation of the convex envelopes in Fig. 10.
< q U be sequences of equidistant points, where q and q are the tangential points defined above. A polyhedral relaxation is given by the following linear inequalities (see Figs. 11a, 12a) :

Fig. 11 Polyhedral relaxations of nonconvex constraints
If N c = 0 , the polyhedral relaxation is defined only by constraints (40) < q U be sequence of equidistant points. In this case, we consider the linear relaxation given by (see Figs. 11b, 12b): If N c = 0 , the corresponding polyhedral relaxation is defined only by constraints (46)-(48).
Analogously, if q L <q < 0 < q U ≤ q and N c > 0 , let q L < q M 1 < ⋯ < q M N c <q be a sequence of equidistant points. We have (see Figs. 11c,12c): When N c = 0 , the polyhedral relaxation is defined by constraints (50)-(52). Assume 0 ≤ q L < q U . Let N c > 0 and consider a sequence of equidistant points q L < q S 1 < q S N c < q U . We define (see Fig. 11d, 12d): If N c = 0 , the polyhedral relaxation is defined only by constraints (54)-(56). Finally, if q L < q U ≤ 0 and N c > 0 , let q L < q M 1 < ⋯ < q M N c < q U be a sequence of equidistant points. A linear relaxation is given by (see Figs. 11e, 12e): When N c = 0 , the polyhedral relaxation is defined only by constraints (58)-(60). The polyhedral relaxations of (26) defined above are illustrated in Figs. 11 and 12. Observe that their strength depends on q L and q U . Smaller ranges for q U − q L lead to tighter linear relaxations. In conclusion, it is possible to write polyhedral relaxations of constraints (2a) as for suitable matrices R t and E t and vectors r t that depend on q t L and q t U .

Appendix 2: An illustrative example
In this appendix, we illustrate the graph decomposition routine outlined in Sect. 3.3 to compute the index set corresponding to links considered by the domain reduction procedure in Algorithm 1. Consider the network ToyNet, whose layout is presented in Fig. 13. This network includes one fixed-head node H 0 , and six unknown-head nodes. Moreover, node V 6 is the only node in the network with degree one. In this case, there is a unique tree, which is composed by links P 6 and P 7 , and nodes V 3 , V 5 , and V 6 . Node V 3 is the root of the tree. Following the notation introduced in Sect. 2, the operation of ToyNet is considered under n l different demand conditions. Let t ∈ {1, … , n l } . Equation (2c) at V 6 implies that q t denotes the demand at node V 6 and time t. Therefore, flow across link P 7 is known a priori, and upper and lower bounds on flow variables corresponding to link P 7 are equal. Furthermore, Eq. (2c) imply that q t P 6 = d V 5 + d V 6 . Hence, also the flow across link P 6 is known a priori. Furthermore, we observe that links P 2 , P 4 and P 5 are connected in series. Equation (2c) yield that As a result, it is possible to select link P 4 as representative and derive bounds on variables q P 2 and q P 5 using (63). In conclusion, the optimization based bound tightening