A Mixed Integer Non-Linear Programming Model for the Optimal Valve Placement within Water Distribution Networks †

: Water leakages represent a crucial aspect in the management of water distribution networks (WDNs). The pressure control, by means of pressure reducing valves, is an economical and viable strategy to contain water leakages and it has been deeply investigated in literature for many years. However, in the presence of high excess pressure and high discharges, thus high energy potential, it may be more convenient to install an energy production device (e.g., a turbine, a pump as turbine) in order not to dissipate the excess pressure, but rather to convert it into energy by means of an electrical generator. Nevertheless, whenever the pressure containment could ensure large water savings but small energy production, the installation of a pressure reducing valve may be more convenient due to the lower purchase and maintenance costs. In literature, many studies have dealt with the optimal location of such valves within WDNs in order to maximize leakage reduction. In this study, the optimization of both the number and location of pressure reducing valves is performed, with the aim of maximizing water savings and minimizing investment cost. Instead of employing external software to compute ﬂow through links and pressure at nodes, the hydraulic resolution of the case study network is coupled with optimization procedure in one single mathematical model. Due to the strong non-linearities modelling of the hydraulic network, as well as the presence of both continuous and integer optimization variables, the resulting mathematical problem is mixed integer non-linear programming (MINLP). A global optimization solver is employed to solve the problem and the achieved results are presented with reference to a literature case study network. According to the results, the proposed optimization ensures an improvement in leakage reduction of 21%, compared to other studies on the same network.


Introduction
Water distribution networks (WDNs) are energy demanding systems affected by very low efficiency [1,2] due to high energy consumption [3,4], as well as large water leakages caused by high pressure values [5,6]. Such a low energy efficiency results in high supply operations costs, as well as in strong environmental impact due to the increase in acid rains and greenhouse emissions [7][8][9]. Hence, the improvement of energy efficiency [10] represents a topic of paramount importance in the management of WDNs. A sustainable growth of such systems may result from strategies to reduce the amount of leaked water, in order to reduce the energy requirements. Despite their effectiveness, the rehabilitation or replacement of damaged pipes are very expensive strategies [11]. A very efficient alternative to contain water leakages is represented by pressure control strategies [12], and installing regulation valves within the network in order to reduce the pressure and, thus, contain the amount of leaked water. In the recent literature, the replacement of valves with energy production devices (EPDs), such as turbines, micro-turbines [13], and pumps as turbines (PATs) [14,15] has been deeply investigated. EPDs, indeed, allow for both pressure 2 of 8 containment and energy production, further improving the energy efficiency of WDNs [16]. Nevertheless, the feasibility of such devices strongly depends on the amount of energy that is recoverable. In the presence of small energy potential (i.e., small head drops and/or small flow rates), energy production may not cover the high purchase and installation costs, thus the installation of pressure valves may represent a more viable solution to save water.
The use of optimization tools for pressure control in WDNs has been significantly investigated in the past and in recent times [17][18][19][20]. The optimal location, as well as the optimal setting of such devices, is crucial for saving the number of leakages and improving the energy efficiency of water systems. In this study, the optimal location of pressure reducing valves (PRVs) is investigated with reference to a literature synthetic network, with the aim of maximizing leakage reduction and minimizing investment costs. The optimization procedure is coupled with the hydraulic resolution of the network in one single mathematical model. A global optimization solver is employed to perform the optimization and the promising results are compared with other literature studies on the same network, in order to show the effectiveness of the proposed optimization.

Optimization Methods: Heuristic and Deterministic Approach
Optimization methods can be classified as heuristic and deterministic approaches. Heuristic approaches are based on computational procedures that search for an optimal solution by iteratively attempting to improve a candidate solution according to a given measure of quality. Despite being very robust, heuristic approaches do not provide any information about the quality of the found solution [21]. Conversely, deterministic approaches rely on the analytical properties of the problem to generate a sequence of points converging to a global optimum or an approximately global optimum. Deterministic approaches include: linear programming (LP), non-linear programming (NLP), mixed-integer linear programming (MILP) and mixed-integer non-linear programming (MINLP) [22]. With reference to the latter problems, it is worth specifying that most MINLP solvers achieve global optima only in convex problems, whereas the solvers managing to find the real optima also in non-convex problems are the global optimization solvers. In this study, the optimization procedure is performed by SCIP (Solving Constraint Integer Programs) solver [23], which is a global optimization solver implementing a spatial branch and bound and various heuristics. It is also worth clarifying that the main reason why a global optimization solver is employed is not the achievement of the absolute optimum (which could require extremely long computational time in hard problems), but rather the fact that this kind of solver provide useful information about the quality of the found solution. For the sake of clarification, with reference to SCIP, during the computation it provides the user with the relative gap between the best solution found and the bound of the problem, which the solver obtains by a relaxation of the problem (Figure 1). During the computation performed by SCIP, the solver reduces the upper bound (or increases the lower bound in case of a minimization problem) and improves the solution, thus the relative gap is progressively reduced until the gap becomes null, which means that the global optimum is achieved.

The Optimization Problem
The analyzed case study is a literature synthetic network [24], which has been considered by many authors in literature. The number of l links and n nodes is equal to 37 and 25, respectively. The network is also fed by three reservoirs, whose level has been assumed as constant. The network layout is presented in Figure 2. The presence of a valve within a link k (k = 1 . . . l) of the network is modeled as a binary variable (I V k ), which is equal to one if the device is installed, zero otherwise. Besides the location, a further variable of the optimization problem is represend by the head-loss within the installed valves, that is, H V k . Since the proposed optimization does not rely on the use of any external softwares to solve the hydraulic network, further variables of the optimization are represented by the discharge flowing through the k-th link (Q k ) and the head at the i-th node (i = 1 . . . n) (H i ). If the node is a reservoir, since the head is known, the variable is the flow rate supplying the network (q r ) at each time interval.
It is worth underlining that, when a deterministic optimization approach is used, the use of non-differentiable functions (e.g absolute value, if function, etc.) should be avoided since these increase the computational complexity of the optimization. As a matter of fact, there exist global optimization solvers that are not even able to handle this kind of functions. Thus, a further binary variable (ζ k ) has been introduced, splitting the flow through links (Q k ) into its positive and negative components [25], depending on the flow direction: where Q max represents the bound of the variable. This formulation allows to avoid the use of non-differentiable functions, which are needed to model the flow direction within the links in the continuity equation at the nodes. The same formulation has been extended to the head drop within the valves, resulting in (H V+ k , H V− k ).

The Mathematical Model
The complete mathematical model describing the optimization problem is presented in Equation (3). Since the variables are both binary and continuous, as well as the con-Environ. Sci. Proc. 2022, 21, 2 4 of 8 straints are both linear and non-linear, the proposed model is mixed integer non-linear programming (MINLP). constraints are both linear and non-linear, the proposed model is mixed integer non-linear programming (MINLP).
With reference to the model in Equation (3), as objective function the net present value (NPV) of the investment has been considered in order to account for both the outflow and inflow cash due to the investment and saved water, respectively. With regard to the investment cost for valve installation, the literature cost model developed by Garcia et al. [26] for PRVs was selected, according to which the total cost of the valves is dependent on pipe diameter, as shown in Figure 3. Moreover, with regard to the water saving (i.e., W in NPV formulation), it has been evaluated as difference between the amount of water that is leaked before performing the pressure control and the leaked volume once the valves are installed within the network.
The leaked flowrate (f p ) has been evaluated according to the formulation proposed by Araujo et al. (2006) [27], depending on the pressure p at the i-th node. Such a leaked flow rate is also taken into account within the continuity equation at each i-th node. Moreover, (H , H , H ) ∈ ℝ, (q , q ) ∈ ℝ, ⩝θ , θ = 1. . n : θ < θ ⩝ i, j = 1. . n, ⩝ k = 1. . n With reference to the model in Equation (3), as objective function the net present value (NPV) of the investment has been considered in order to account for both the outflow and inflow cash due to the investment and saved water, respectively. With regard to the investment cost for valve installation, the literature cost model developed by Garcia et al. [26] for PRVs was selected, according to which the total cost of the valves is dependent on pipe diameter, as shown in Figure 3. Moreover, with regard to the water saving (i.e., W in NPV formulation), it has been evaluated as difference between the amount of water that is leaked before performing the pressure control and the leaked volume once the valves are installed within the network. The leaked flowrate (f p ) has been evaluated according to the formulation proposed by Araujo et al. (2006) [27], depending on the pressure p at the i-th node. Such a leaked flow rate is also taken into account within the continuity equation at each i-th node. Moreover, (H , H , H ) ∈ ℝ, (q , q ) ∈ ℝ, ⩝θ , θ = 1. . n : θ < θ ⩝ i, j = 1. . n, ⩝ k = 1. . n With reference to the model in Equation (3), as objective function the net present value (NPV) of the investment has been considered in order to account for both the outflow and inflow cash due to the investment and saved water, respectively. With regard to the investment cost for valve installation, the literature cost model developed by Garcia et al. [26] for PRVs was selected, according to which the total cost of the valves is dependent on pipe diameter, as shown in Figure 3. Moreover, with regard to the water saving (i.e., W in NPV formulation), it has been evaluated as difference between the amount of water that is leaked before performing the pressure control and the leaked volume once the valves are installed within the network. The leaked flowrate (f p ) has been evaluated according to the formulation proposed by Araujo et al. (2006) [27], depending on the pressure p at the i-th node. Such a leaked flow rate is also taken into account within the continuity equation at each i-th node. Moreover, With reference to the model in Equation (3), as objective function the net present value (NPV) of the investment has been considered in order to account for both the outflow and inflow cash due to the investment and saved water, respectively. With regard to the investment cost for valve installation, the literature cost model developed by Garcia et al. [26] for PRVs was selected, according to which the total cost of the valves is dependent on pipe diameter, as shown in Figure 3. Moreover, with regard to the water saving (i.e., W v y in NPV formulation), it has been evaluated as difference between the amount of water that is leaked before performing the pressure control and the leaked volume once the valves are installed within the network. The leaked flowrate (f i p β i ) has been evaluated according to the formulation proposed by Araujo et al. (2006) [27], depending on the pressure p at the i-th node. Such a leaked flow rate is also taken into account within the continuity equation at each i-th node. Moreover, besides the leaked flow rate and discharges flowing through the links approaching the i-th Environ. Sci. Proc. 2022, 21, 2 5 of 8 node, the continuity equation also accounts for the water demand of the node, expressed as follows: where q d i is the daily average demand of i-th node and c d (t) is the demand coefficient at time t, varying according to the pattern used in [17].
Moreover, with regard to the momentum balance constraint, r k L k is the head loss within the k-th link with length L k , and it is evaluated according to Hazen-Williams formula: With reference to the model in Equation (3), the continuity equation at the nodes and the momentum balance equation within the links have been expressed as inequalities, since two small feasibility tolerances (tol Q i and tol H k ) have been introduced in order to improve the convergence of the solver. The tolerance formulations are following presented: where ε Q is equal to 0.01, and J i is the set of the j nodes connected to the i-th node by a single pipe. With regard to the expression of tol H k , ρ k has been already defined in Equation (5) and ε H has been assumed equal to 100. However, both the resulting tolerances represent very small errors, which can be considered as acceptable as the errors resulting from the Hazen-Williams resistance formula or water demand estimation.
As shown in the mathematical model, a further constraint is represented by the pressure at nodes (i.e., H i − z i ), which has been bounded between a minimum and a maximum value (i.e., p min γ and p max γ , respectively, where γ is the specific weight of water), set as 25 m and 100 m, respectively. Moreover, in order to avoid the installation of valves where the flow reverts, a set of constraints [17] has been defined, forcing the variable I V k to be equal to zero when the variable ζ k changes at least once during the day, among all possible combination of time intervals θ 1 (θ 1 = 1..n d ), θ 2 (θ 2 = 1..n d ), with θ 1 ≤ θ 2 : Finally, a maximum value for the variable head-loss has been defined (H k max ) in order to bound the variable when the valve is installed (i.e., I V k is equal to 1), and force it to be equal 0 otherwise. However, such a constraint is not intended for its proper hydraulic meaning, but rather has been defined for reducing the research space and enhancing the convergence of the problem.
It is worth underlining that, despite the fact that in the model formulation the dependence on the time interval is omitted for the sake of brevity notation, all the variables and constraints are computed along the time. The total number of variables amounts to 1545 (i.e., 296 binary and 1249 continuous variables), whereas the total number of mathematical constraints is equal to 4051.
Finally, the optimization code has been developed by A Mathematical Programming Language (AMPL) [28] and the procedure has been implemented on an Intel @ Xeon(R) CPU E5-2620 v4 @ 2.10 GHz × 16 with 64 GB RAM.

Results
The optimization results are presented in Table 1. According to the found solution, the NPV amounts to EUR 709,186. The solver selects nine valves requiring an investment cost equal to EUR 72,361 and ensuring a water saving of 924 m 3 per day. With regard to the solution quality, the relative gap provided by the solver is equal to 3%, which means that the solution is approximately the global optimum. It is worth clarifying that the solver struggles more to decrease the bound of the problem than to improve the incumbent solution. Therefore, even when the global optimum is found, the gap could be not perfectly zero. However, in MINLPs affected by hard computational complexity, the achievement of the 0% gap may require even an infinite time. Table 2 shows a comparison with other literature studies on the same network. Nicolini and Zovatto (2009) [19] performed the optimization by using multi-objective genetic algorithms. In particular, the first criterion was represented by the minimization of the number of valves, and the second was the minimization of the total leakage in the system. Dai and Li (2014) [20], instead, minimized the total excessive pressure, reformulating the MINLP problem to a mathematical program with complementarity constraints, which can be efficiently solved by available NLP algorithms. In both studies, the optimization was performed by considering three load conditions, corresponding to demand multiplier values (c d ) as 0.6, 1.0 and 1.4. According to Table 2, compared to the literature studies, the proposed optimization selected a higher number of valves (i.e., 9 against 5). Regarding the location, it essentially depends on the employed cost model. Indeed, in the proposed work, the solver selects the links with small diameters in order to minimize the valve cost according to Figure 3. Moreover, with regard to the average leakage (i.e., the total leaked flow rate averaged on the three load conditions), in this study it is significantly smaller than the amount obtained in both the literature studies and this results from the installation of a higher number of devices. However, given the cost model of Figure 3 and assuming the NPV as the objective function, the result achieved by the proposed optimization is a near global optimum. Therefore, beyond the promising results, the effectiveness of the proposed optimization lies in the evaluation of the quality of the found solution, which is not possible when heuristic approaches are used.

Conclusions
In this study, the optimal location of pressure reducing valves is investigated within a literature synthetic network widely used in literature. A deterministic approach has been used with the aim of finding the global optimum (or near global optimum) of the optimization problem. A mathematical model has been developed, consisting of both the optimization procedure and hydraulic equations to compute the flow through links and pressure at the nodes. To perform the optimization, SCIP solver has been used, which potentially finds the global optimum in MINLP problems, no matter their convexity. According to the results, a very good quality solution has been obtained, presenting a relative gap value with the bound of the problem as 3%. The solution consists of an NPV as EUR 709,186 and the resulting water saving amounts to 924 m 3 per day. The solution has been compared with other literature studies in terms of position and average leakage. According to the comparison, the proposed optimization selects a higher number of valves (i.e., 9 against 5), ensuring a more significant leakage reduction. However, the real strength of the proposed optimization lies in the possibility for the user to evaluate the effectiveness and percentage of improvement of the found optima, which is not allowed when heuristic approaches are employed to perform the optimization.