A MINLP model for optimal localization of pumps as turbines in water distribution systems considering power generation constraints

: Pressure reducing valves (PRVs) are commonly used for pressure control in water distribution systems (WDSs) by means of dissipating the pressure excess. The use of pumps as turbines (PATs) is an alternative and more favorable system since they not only control the system pressure to decrease water leakage, but also utilize the pressure excess to generate electrical energy. The optimal localization of PATs can be casted into a mixed-integer nonlinear program (MINLP) where binary variables are used to represent the presence of PATs on links. Most of the available MINLP models for optimal PAT localization adopted the optimization approaches for PRV localization without considering the bound constraints on ﬂow rates and heads of PATs. As a result, such an optimization model may make PATs delivering a non-desired output. In this paper, we propose a new MINLP model for optimal PAT localization. Instead of using a constraint on the maximum number of PATs to be placed in a WDS, new constraints relating to the minimum power generated by PAT are introduced to ﬁnd links having adequate ﬂows and head drops for placing PATs. Moreover, constraints are used to restrict ﬂows and heads of PATs to their feasible operating range, so that the problem can be e ﬃ ciently solved. The proposed MINLP model is applied to the optimal localization of PATs for a WDS benchmark and a real-world WDS in Vietnam. The results demonstrate that the new MINLP model can e ﬃ ciently identify optimal locations for PAT placement where the speciﬁed working range and minimum power generated by the PATs are ensured.


Introduction
Water loss occurs in all water distribution systems (WDSs). The main cause of water loss is due to inefficient pressure management. Many strategies have been applied to decrease the water loss, such as rehabilitation/repair of WDSs including replacement of pressure release valves, early detection of water leakage/break points, and pressure control [1,2]. Due to the fact that water leakage is proportional to pressures at nodes, the control of the pressure in WDSs is considered as one of the most efficient and cost-effective approaches to decrease water leakage and the probability of newly created leakage points [3]. Pressure reducing valves (PRVs) are commonly installed in WDSs to control the system pressure for water leakage reduction by dissipating the potential energy of the water in WDSs [4]. In the literature, there are many optimization models and solution approaches to solve the problem for PRV localization. Eck and Mesvissen [5] developed a mixed-integer nonlinear program (MINLP) model for identifying optimal locations for PRV placements in WDSs. Dai and Li [6] reformulated the MINLP model for optimal PRV localization in [4] into a mathematical program with complementarity constraints (MPCC), allowing standard nonlinear programing (NLP) techniques applied to solve the MPCC efficiently. Ajauro et al. [2] employed a genetic algorithm (GA) to address the optimal locations for placing PRVs. NSGA-II and scatter-search meta-heuristic algorithms were also applied to solve PRV locating problems [7,8]. Pecci et al. [9] and Pecci et al. [10,11] proposed a MINLP model and an efficient solution approach to optimal PRV localization by significantly reducing the number of nonlinear constraints in the MINLP. Recently, Cao et al. [12] proposed a method to simultaneously place pressure sensors and localizing PRVs for the purpose of control system design for on-line pressure control.
Although PRVs are commonly installed in WDSs to control water leakage, they dissipate the head loss excess without any profit. Fontana et al. [13] proposed a method to control the water leakage and simultaneously recover electrical energy by means of installing hydropower stations in WDSs. For the task of energy recovery, the application of pumps as turbines (PATs) is considered as an emerging hydropower technology for managing WDSs by coupling pressure control with hydropower generation [14,15]. To realize energy recovery in this way, the first task is to determine optimal locations for PAT placement in WDSs, which can be casted into a mixed-integer nonlinear program (MINLP) [16,17]. In addition, the mechanism of hydraulic regulation (HR) or electrical regulation (ER) can be used to maximize the energy production by the application of the variable operating strategy (VOS), as studied in [18,19]. The HR moves the available pressures and flows to the ones in the PAT characteristics by regulating the flow control valve and pressure reducing valve in such a way that the PAT system efficiency is maximized. In contrast, ER adjusts the speed of the generator in the PAT system to fit their characteristics to the available flows and head drops. Previous studies on employing PATs in WDSs can be roughly classified into economic analysis and solution approaches to optimal localization.

Studies on Economic Analysis by Using PATs for Energy Recovery
For the potential analysis of installing PATs for energy recovery and pressure management, Puleo et al. [20] investigated the application of centrifugal PATs in a real WDS for potential energy recovery with the presence of private tanks and intermittent service. The results show that the energy production may be low and discontinuous, and hence, the efficiency of the PAT installation in WDS should be studied. Carravetta et al. [21] evaluated the beneficial cost for installing PATs in WDSs by proposing an efficient control scheme for bypass and control valves. Fecarotta et al. [22] developed a so-called business plan model for the design of the hydropower plant. Based on that, the economic benefit of valve replacement by PATs can be evaluated.
De Marchis et al. [23] analyzed the energy recovery in the Palermo network in Italy. The results indicate that the application of PATs can lead to a very attractive economic benefit. Samora et al. [24] proposed a method to quantify the potential for hydropower, especially for micro-turbines with a five-blade tubular propeller. The method was evaluated for the WDS in the city of Fribourg, Switzerland, demonstrating that about 10% the energy potential can be recovered. Parra et al. [25] proposed a system composed of a PRV and a PAT in combination with intelligent pressure management to reduce water loss and recover the energy. The operational combination of a PAT and a PRV has shown that the water saving and electrical energy recovery are increased and that the PAT-PRV-system is suitable for WDSs with high difference elevations, high operational pressures, and high demand variability.

Studies on Optimal Localization of PATs by Stochastic Search Approaches
For optimal localization of PATs to maximize energy recovery and leakage reduction, Giugni et al. [26] proposed a new formulation of the objective function for optimal PAT localization. Instead of minimizing the excessive pressure, the electric power generated by PATs is maximized. The optimization problem was solved by using a genetic algorithm (GA) and the results show that by using this objective function, although the amount of water leakage is slightly larger than that resulted by using the minimization of excessive pressure, the energy recovery is much higher.
Jafari et al. [27] proposed two steps to control pressure and recover energy. In the first step, optimal placements and settings of PRVs were determined to minimize the excessive pressure by using GA. In the second step, PRVs with a high head loss and adequate flows were replaced by PATs. The EPANET toolkit [28] was applied to carry out simulation of the WDS for the GA implementation. The solution approach was applied to both pressure control of a real WDS in Iran and energy recovery, where it is shown that by replacing PRVs by PATs, an excessive pressure decrease can be achieved in combination with energy recovery.
Lima et al. [29,30] presented a method for PAT selection and localization in WDSs based on maximizing the cost of energy recovery and the saving cost of water leakage reduction where the PATs speeds can be varied to improve the efficiency. Constraints were used to maintain the pressures at nodes within their bounds. Particle swarm optimization (PSO) was applied to solve the optimization problem, while rotational speeds of PATs were calculated according to the nearest curve available in a set of complete characteristic curves of the pumps in the Suter plan. Tricario et al. [31] optimized the energy recovery by PATs for improving the WDS operation by formulating a multi-objective optimization problem where the decision variables are pump scheduling, locations, and types of PATs, as well as initial water tank levels. The terms in the objective function were pumping energy costs, pressure surplus, and PAT recovery energy cost. The trade-off results allow users to select appropriate pump scheduling so that water leakage in the WDS is controlled while the operating cost can be significantly reduced with some pumps working in the turbine mode.

Studies on Optimal Localization of PATs by MINLP Approaches
Corcoran et al. [16] optimized the operation of WDSs for the combination of hydropower energy recovery and leakage reduction. The optimization problem was formulated and solved for determining optimal PAT locations to maximize the electric power recovery. The solutions from various optimization models such as MINLP, nonlinear program (NLP), and GA were compared. The comparisons revealed that the solution by the MINLP model was the most suitable for determining optimal PAT locations. Fecarotta and McNabola [17] proposed another MINLP model for the optimal localization of PATs in WDSs, in which the profitable cost obtained with energy recovery and leakage reduction was maximized. Although the optimal solution balances the two objectives, the formulation of the MINLP has some weaknesses. First, the flow through the PATs can be reverse according to the change of the demand during the day, which is not possible in practice. Second, there is no constraint relating to the minimum powers that PATs will generate, i.e., power generated by PATs can be very low.
In summary, most previous works relating to pressure regulation and energy recovery in WDSs have focused on evaluation of economic aspects with replacement of PRVs by PATs or optimal localization and operation. In addition, most of the optimization models used for the optimization of PAT localization were adopted from the existing models used for optimal PRV localization, e.g., the MINLP model in [16] used in [5] or the optimization based EPANET simulation model used in [26]. The differences between these models only come from the definition of the objective function in the optimization model, i.e., maximization of the power production used in [16,26], or profit maximization from energy recovery and cost minimization of water lost in [17] for the same optimization problem. From the practical operation point of view, with PRVs, there is no requirement on flows and head drops (i.e., a PRV can be operated in a closed mode). In contrast, with PATs, since the flows and head drops are related to their generated power, their operation must be bounded to the specified feasible ranges (i.e., a PAT cannot be operated in any values of flows and head drops). Such a boundary value is crucial for PATs avoiding low energy production and discontinuous operations [32,33]. In addition, in the design of a WDS with PATs, the predominant aim is to achieve the practically relevant system pressure for water supply. However, to the best of our knowledge, there is only one study [17] concerning these issues, where only the lower bounds of head drops of PATs were considered. Therefore, in this paper, as considered in [26], we address the problem of selecting optimal PAT locations, aiming to maximize the electric power production while regulating the system pressure for water leakage reduction. For this purpose, a new MINLP model is proposed with constraints relating to the specified minimum powers generated by PATs and constraints on the bounds of flows and heads across the PATs. In this way, the required minimum powers generated by PATs are ensured and the PATs to be placed will be operated in the practically allowable region. It means that only links that have adequate flows and head drops will be chosen for placing PATs. In addition, the optimization model is suitable for PAT operations in practice where reverse flows will be prevented. A benchmark and a real-world WDS are used to demonstrate the usefulness of the proposed model.
The remainder of this paper is organized as follows. In Section 2, the new MINLP model for PAT localization is formulated, where new constraints relating to the power generation of PATs and bounds of flows and heads are proposed. Two case studies are carried out in Section 3 and conclusions are given in Section 4.

Problem Formulation for Optimal PAT Localization
The aim of this study is to develop an optimization model for identifying optimal locations for placing PATs with which the total energy generation is maximized and at the same time the excessive pressure reduced. Since the produced energy depends on flows and head drops across PATs, optimal positions of PATs will achieve a maximum energy recovery. In addition, as the PATs absorb much more excessive pressure, it will lead to a high reduction of the water leakage in the WDS.
Since PATs can be placed on a link in either forward or backward flows, we define two binary variables, v i,j and v j,i , for each link between node i and node j to indicate whether a PAT is placed on the link in forward (from i to j) or backward (from j to i) flows, respectively. In particular, if v i,j = 1, the PAT will be placed on link ij in forward flows and otherwise if v j,i = 1, the PAT will be placed on link ij in backward flows. Since in each link it is possible to place only one PAT in one direction of flow, v i,j and v j,i cannot be accepted for the value of 1 simultaneously; this will be ensured by a constraint (i.e., Equation (13)). In addition, we introduce a continuous variable θ i,j,k for each link describing the head drop across the PAT.

Objective Function
We consider a WDS with NP links, NJ nodes, NR reservoirs, and NL demand patterns. The objective function, defined in the optimization problem, is to maximize the total net power generated by the PATs installed in the WDS. The power generated by PATs depends on the flow Q i,j,k and head drop across the PAT (θ i,j,k ). Therefore, the objective function is defined as: where Q i,j,k is the flow through the PATs; γ = 9806 N/m 3 is the specific weight of water; and η i,j,k is the average efficiency of the PATs at time step k, it is taken a value of 0.65 as used in [17]; k = 1, . . . , NL is the time step. In this study, by assuming a constant efficiency, as in [17,26], our optimization model provides the average energy generated by the PATs, i.e., the results present a reference value for the energy recovery. It is noted that, for a more realistic implementation, the efficiency of the PATs should be considered by integrating their characteristic curves into the optimization model, i.e., the PATs' efficiency (η i,j,k ) varies with the flows according to their efficiency curves. The effect of variations in PATs' efficiency with respect to variations in flow and pressure, due to a diurnal demand pattern, can be accessed when PAT types are chosen [16]. In addition, the method of variable operating strategy (VOS) proposed in [18] can be applied to match the flows and head drops across the PATs with their characteristic curves. These will be considered in an extended optimization model in our future research.

Constraints
For a WDS with its layout and demand patterns given a priori, the optimal decision for PAT localization is to determine optimal locations for placing the PATs at which the maximum amount of power can be generated. In addition, operational constraints for efficient operations of the WDS and PATs must also be met, such as maintaining adequate services with adequate water pressures and flows. In addition, the PATs should be operated in their specified working ranges [32,33]. Therefore, the optimization problem for PAT localization is subject to the following constraints.
(1) Continuity flow The continuity equation at node i for scenario k: where Q j,i,k and d i,k are the flow through link ji and demand at node i at time step k, respectively; the leakage amount l i,k associated to node i is calculated by [2,6]: where L t,i is the total length of links connected to node i, calculated by: where E i is the elevation of the node; p i,k = H i,k − E i and H i,k is the static pressure and total head at node i, respectively; β is the leakage exponent, C L is the discharge coefficient of the orifice and has the unit of l/ s · m 1+β [2,17]; L i,j is the length of link ij. The units of elevation, pressure, head, and link length are meters.
(2) Conservation of energy The hydraulic model equation describing whether a PAT is installed on the link between node i and node j is written as [9][10][11]: where ∆H i,j,k is the head loss across the link; H i,k and H j,k are the nodal heads at node i and j at time step k, respectively; θ i,j,k is the head drop across the PAT. When θ i,j,k = 0, no PAT is placed on the link, and Equation (5) 1.852 (6) or by the Darcy-Weisbach equation: where L i,j , D i,j , and C i,j in (6) and (7) are the length, diameter, and Hazen-William coefficient of link ij, respectively; f i,j in (7) is the pipe friction factor; g is the acceleration of gravity.
(3) Bound constraints for flows and head drop across PATs According to different types of PATs, it is necessary to ensure that their head drops must be bounded [17,32,33], namely θ i,j,k ≥ H min T . This requirement can be accomplished by the following constraints: where H min T is the minimum head required for placing a PAT, the value of which can be properly chosen according to the type of PATs or expected lower bound of the heads for installing the PATs. Similarly, the flows through the PATs should be in their feasible range [32,33], i.e., Q Tmin for backward flow. We use the following constraints to ensure these requirements: where Q Tmin i,j is the minimum flow of the PATs or specified range of the flows for installing the PATs. The value of Q max i,j is the maximum allowable velocity through the pipes. The constraints in (10) and (11) exactly satisfy the bounds of flows. In fact, when v i,j = 1 and v j,i = 0(i.e., a PAT is placed on link ij with direction from i to j), the constraint in (10) becomes when v i,j = 0 and v j,i = 1 (i.e., a PAT is placed on link ij with direction from j to i). When v i,j = 0 and v j,i = 0, the constraints in (10) and (11) (4) Constraints on the minimum power generation Instead of fixing the number of PATs to be placed in the WDS, our aim is to find links for placing PATs from which a specified minimum amount of the generated power will be ensured, e.g., larger than P min given a priori. The idea of ensuring the minimum power generated by PATs comes from [17] due to the generator specification. However, the authors in [17] used an empirical method (i.e., after the MINLP algorithm identifies the optimal locations of PATs, only the PATs generating powers larger than P min are chosen for the installation), instead of using constraints to eliminate links providing power lower than P min . Similarly, the authors in [16] used a nonlinear programming method, instead of MINLP, to determine optimal locations of PATs where all links in WDS are considered as candidates, but only links providing power higher than P min will be chosen. In our study, we introduce this constraint on the minimum power generated by PATs explicitly in the MINLP framework. The choice of the minimum power value P min depends on the type of PATs to be installed and will prevent the PATs from operating at low power [32]. Such a constraint is expressed as: where s i,j,k ≥ 0 is a tolerant variable; its upper bound can be set to a small value, e.g., 0.01.
The introduction of such a variable allows links capable of delivering the minimum power with a minor violation.
(5) Constraints on binary variables associated to each link To ensure that a PAT can only be placed in either forward or backward flows, it should be: From the PAT operation point of view, both the specified minimum energy and the flow as well as head constraints must be satisfied. From the mathematical solution point of view, all constraints are independent and there is no priority between them, i.e., the solution of the optimization problem satisfies all the constraints.
In addition, the head drops across the PATs are defined as optimization variables in our MINLP model, and an equation (i.e., Equation (12)) is used to explicitly satisfy the minimum power generation by the PATs. This is not the case in the MINLP model in [5] or the optimization-based simulation model used in [25] where the head drop of PAT is implicitly included in the deviation of the heads between the upstream and downstream nodes. Comparing with the MINLP model in [17], our model will be more efficient because of preventing reverse flows and ensuring the minimum power generated by the PATs. To the best of our knowledge, such a formulation of the MINLP model has not been applied to solve the optimal PAT localization problem.
In the next section, we will evaluate the performances of our MINLP model with a variety of scenarios on the specified minimum power. A MINLP solver, BONMIN [35] in GAMs [36], is used to solve the corresponding problems. The branch-and-bound algorithm is implemented in BONMIN to solve the MINLP problem. This algorithm is based on solving a continuous nonlinear program at each node of the search tree and branching on the binary variables [35].

Case Studies
We consider two case studies where the first one is the WDS benchmark studied in [26], and the second one is a real WDS in a City in Vietnam for optimal PAT localization. The advantage of our MINLP model is that constraints to satisfy the lower bounded power generated by PATs and constraints on flows and head drops of the PATs are introduced in the problem formulation. To demonstrate the effects of these constraints, in each case study, we apply our approach to determine the number of PATs and their locations with respect to several scenarios with different lower bound values specified to the PATs. Beside considering the benchmark WDS studied in [26], which is a small scale WDS with a small number of links and nodes, a real large-scale WDS in a city in Vietnam is also taken to demonstrate the capability and efficiency of our approach for complicated WDSs.

Case Study 1
We first consider a benchmark WDS, as shown in Figure 1, comprising of 37 links, 22 nodes, and three reservoirs, which was studied in [1,2] for optimal pressure management. The data of links and nodes are given in [1,2], while the discharge coefficient C L is 10 −5 and the leakage exponent parameter β is 1.18 [2], respectively. The head losses across links are calculated by using the Hazen-Williams equations (Equation (6)). The minimum allowable pressure at nodes is 25.0 m for the water supply [17,26]. The 24 demand pattern factors representing the hourly changes of demands at nodes and 24 reservoir water levels taken from [2] are given in Table 1 and Figure 2a, respectively. In addition, for the purpose of comparison with the results reported in [26], the normal reservoir head values used in [17,26] are also considered, i.e., the water levels of reservoirs 23, 24, and 25 are 55.6 m, 55.5 m, and 56.0 m, respectively.   = 600 L/s). The average efficiency of the PATs is taken as 0.65 [17].
To demonstrate the efficacy of approach, we consider a variety of scenarios by restricting the minimum power generated by the PATs, i.e., P min is defined as 0.25 kW, 0.75 kW, and 1.0 kW, respectively. The results of optimal locations for placing PATs and computation time required for solving the optimization problem are given in Table 2. According to different minimum power values, we found different numbers of links suitable for placing PATs. For the first scenario (i.e., P min = 0.25 kW), four links suitable for placing PATs are determined, which are 23-1, 13-12, 25-16, and 24-10, and it took 3542.20 s to solve the problem. The hourly powers generated by PATs are illustrated in Figure 2b, showing that the PATs indeed generate powers higher than 0.25 kW as expected. The total electrical energy recovery per day will be 221.19 kWh. As shown in Figure 3a,b, the flows and head drop across the PATs are higher than 10 L/s and 4.0 m, respectively, thus, satisfying the specified bounds for the PATs. It is also seen that, when a high demand appears, i.e., from time 9:00 to 18:00, the flows through the PATs also increase, while the head drops across the PATs decrease. This is because that the reservoir water levels change slowly while the demand patterns vary highly, i.e., during the time from 4:00 to 10:00 as seen in Table 1. Moreover, among the four PATs determined, the PAT on link 25-16 produces the lowest amount of electric power, while the PAT on link 13-12 generates the highest.
In scenario 2, P min = 0.75 kW, the optimal solution found three links, 13-12, 24-10, and 23-1 for placing PATs and the computation time was 1537.86 s. As seen in Figure 4a, the powers generated by the PATs in each hour are higher than 0.75 kW. More importantly, once again we observe that the PAT on link 13-12 produces the highest amount of electric power, while the lowest will be generated by the PAT located on link 23-1. To assess the effect of the variation of the reservoir water levels, we computed the electrical energy recovery in the case where the reservoir water levels change during 24 h as in [2] and the case where they are fixed to normal values as in [26], respectively. It is seen from Figure 4b that there is a slight difference in the energy recovery and leakage reduction.  Interestingly, our approach found the same optimal locations for three PATs as by Giugni et al. [26]. By comparison, the potential electrical energy reported in [26] needs to be multiplied by a factor of 0.65 representing the turbine efficiency, which was also used in [17]. The comparison on the electrical energy recovery is given in Table 3. Moreover, the comparison on the hourly generated energy is shown in Figure 4b. In addition, in the case of constant reservoir water levels (i.e., fixed to normal values), the generated energy reported in [26] is a bit higher than that from our approach. The reason is that the electrical energy in [26] was calculated based on the head loss of the link where the PAT is placed (e.g., ∆H i,j,k ), while our approach uses the head drop across the PAT, i.e., θ i,j,k for calculating the electrical energy. For this reason, with the same flows, the energy calculated in [26] will be larger, but this calculation is in fact not accurate. Additionally, it can be seen in the Figure 4b that the energy generated by the PATs will be high when flows increase, and the reservoir water levels are at moderate values. For example, at 9:00 am, the energy attains the maximum value.  Figure 5a it is seen that the flows though PATs are larger than 10 L/s, i.e., satisfying the constraint of allowable flows. Moreover, it shows that the PAT on link 13-12 is suitable for a type of PAT working in a range from 40 L/s to 130 L/s, while the PAT on link 24-10 has a working range from 10 L/s to 30 L/s. The PAT on link 23-1 generates the least power with a working range of flows from 10 L/s to 20 L/s. The head drops across PATs on link 13-12, 24-10, and 23-1 is shown in Figure 5b, showing that the head drops are in a range from 9.0 m to 14.5 m. From these results, it can be concluded that the constraints on flows and head drops introduced in this study make the MINLP model more practical, i.e., for a given type of PAT, it is necessary to restrict the flows and head drops across the PAT to the specified operating range despite of varied demand patterns. In this way, feasible operation of PAT is ensured, and low energy efficiency can be avoided. In other words, for a given WDS, with our approach, it is possible to select links providing the suitable range of flows and head drops for types of PATs.
In scenario 3, P min = 1.5 kW, we found two links, 13-12 and 24-10, for placing PATs, and the hourly powers generated are shown in Figure 6a. The optimal locations of the 2 PATs are the same as those from Giugni et al. [26] in the case of two PATs. The comparison of the electrical energy recovery is given in Table 3 and the hourly energy profiles by both approaches are shown in Figure 6b. It is seen that both profiles are nearly the same. The flows through the PATs shown in Figure 7a are higher than 10 L/s, while the head drops across all PATs take values from 9.0 m to 15.0 m. Based on that, a type of PAT can be chosen, satisfying such operating range with high efficiency of the generators [32,33]. Concerning the water leakage reduction, we compare the leakage flows reported in [26] and those by our MINLP model in Figure 8 for the case of constant heads of reservoirs. It is seen that using our approach, the resulting water leakage reduction is a bit higher than that in [26]. With two PATs, our approach leads to an average reduction of 6.70 L/s while the approach in [26] to 6.51 L/s. Similarly, with three PATs, an average of 7.88 L/s from our model and 7.72 L/s from [26] can be observed. In addition, as shown in Table 3, for four PATs, as compared with the case of three PATs, more energy recovery can be achieved but with a bit lower water leakage reduction. From these results using our MINLP model, the number of PATs to be placed in the WDS depends on the specified minimum power value, P min . As a result, if the type of PATs, the expected operating ranges, and the minimum power requirement are given, our approach can identify links satisfying these conditions for placing PATs. In this case study, we apply our MINLP model to identify optimal locations of PATs in different scenarios on the minimum power setting values. Our contribution is on ensuring the minimum power generated by PATs and keeping their flows and head drops inside the specified bounds. For scenarios where the minimum power is set to 0.75 kW and 1.0 kW, the MINLP model finds the same optimal locations for placing PATs as reported in [26] for the case of three PATs, and two PATs, respectively. Unlike the approaches in [17,26], our approach ensures that the powers generated by the PATs be higher than the specified lower bounds. For the case of three PATs, the generated powers are higher than 0.75 kW, while for the case of two PATs, they are higher than 1.5 kW.

Case Study 2: Optimal Pressure Management for a Real-World WDS in a City in Vietnam
In this case study, we apply our MINLP model for identifying optimal locations of PATs in a large-scale WDS in Thainguyen, Vietnam.
The WDS, shown in Figure 9, consists of 116 pipes, 99 nodes, and four reservoirs. The four reservoirs, located at high elevations, supply water for the residents in the city through the pipeline system. Like case study 1, the Hazen-Williams equations (Equation (6)) are used for the head loss calculation. Although it has not been realized yet, pumps as Turbines (PATs) can be considered in the WDS to simultaneously control the system pressure and recover electrical energy. Therefore, we apply our MINLP model to investigate the potential of energy recovery. In this case study, the network structure and pipeline data are from the local water company. The water demand profiles are from the historical data of the network operation as shown in Table 4   Like case study 1, we consider three scenarios on the minimum power values (P min ) generated by PATs, which are 5 kW, 8 kW, and 12 kW, respectively. In practice, these values are determined by the types of PATs available. As an example, with a centrifugal PAT, a double suction PAT, or a parallel centrifugal PAT, the power generation can vary from 1 kW to 100 kW [32]. The minimum head drops across PATs are set to 2.0 m for all scenarios. The BONMIN solver [35] was used to solve the MINLP problem for these scenarios. The results of PAT locations and computation times are given in Table 5. After the optimal locations of PATs are found, we solve the continuous nonlinear optimization problem (i.e., with the fixed locations of the PATs) for the 24 demand patterns to evaluate the benefit of placing PATs by reducing the excessive pressure. The average power (P av ) and average excessive pressure (EX) are calculated as follows: (18) where N n is the number of demand nodes for which p i,k ≥ p L i ; NL is the number of demand patterns; N PAT is the number of PATs. The results of each scenario are given in Table 5. In the first scenario, P min = 5.0 kW, four links are found for placing PATs, i.e., 151-34, 152-29, 3-63, and 49-47. The corresponding profiles of hourly generated powers are shown in Figure 10. All PATs generate powers higher than 5.0 kW, satisfying the minimum power requirement. The produced energy per day is 1958.0 kWh, where the highest energy is generated by the PAT placed on the link 49-47. The head drops and flows through PAT are given in Figure 11a,b. It is seen that flows are in the range from 50 L/s to 650 L/s, while the head drops belong to the range from 2.0 m to 35.0 m. Based on these results, it is possible to choose suitable PATs [32]. The average excessive pressure is 12.77 m, i.e., reducing nearly a half of excessive pressure amount as compared with the case of no PAT where the average excessive pressure is 23.73 m.  In scenario 2, P min = 8.0 kW, solving the problem results in three links suitable for placing PATs. The hourly powers of these PATs are shown in Figure 12. The produced energy is 1585.0 kWh/day reducing 373.0 kWh/day as compared with the case of 4 PATs. This is reasonable because the number of PATs will decrease if we increase the boundary value of power (P min ). The head drops and flows through the PATs are given in Figure 13a,b, respectively. It is seen that the head drops are between 2.0 m and 35.0 m, while the flows are in the range from 70 L/s to over 600 L/s. The head drops across PATs vary according to the changes of demand patterns. They are high in the low demand periods while they are low in the high demand periods. In all cases, the head drops are bounded to be larger than 2.0 m. The average excessive pressure is 12.79 m, thus reducing nearly 10.0 m as compared with the case where no PAT is used. In scenario 3, P min = 12 kW, we found only one link (41-40) suitable for placing PATs. The average power produced is 48.52 kW and the average excessive pressure is 17.35 m. As compared with the case of no PAT, the average excessive pressure reduces about 6.0 m.
In order to see the reduction of excessive pressure, the pressure contours of the WDS for the case of no PAT (left) and the case of 4 PATs (right) corresponding to the demand patterns of 0.36 and 1.34 are shown in the Figures 14 and 15, respectively. It is seen that a significant amount of excessive pressure is absorbed by the PATs, especially in the low demand patterns. The pressure contours are gained based on simulations with EPANET 2 [37]. In addition, the cumulative distributions of the pressures in the WDS for different demand patterns in 24 h are shown in Figure 16 for all three scenarios. When considering the case of four PATs installed, it is seen that about 60% of the nodes has the pressures lower than 30.0 m and nearly 90% of the nodes with pressures less than 40.0 m. Similar observations can be seen in the case of three PATs. In the case of one PAT, the excessive pressure reduction is less than that in the cases of three and four PATs, but as compared with the case of no PAT, much more excessive pressure is absorbed by the PAT. Based on the above results, the relations between the numbers of PATs to be placed in the WDS, the generated power, and average excessive pressure can be determined to support water utilities for making decisions on the investment of installing PATs. In addition, the proposed MINLP model is general and can be applied for any WDSs for optimal PAT localization. The solution of the optimization problem provides potential links for placing PATs, which will generate a certain amount of power. Based on these results, together with the characteristics of PATs and their efficiency curves, the variable operating strategy (VOS), as studied in [18,19], can be applied to operate the PAT systems in such a way that the overall efficiency will be maximized.

Conclusions
This paper proposed a new MINLP model for optimal PAT localization. Constraints relating to the specified minimum power of PATs are introduced so that the formulated MINLP is capable of determining links having adequate flows and head drops for placing PATs, with which the power generated by PATs is ensured in the specified range. This will avoid the problem of low energy generation and discontinuous operation of PATs. Compared with available MINLP models, our MINLP model is more practical in that the flows and head drops across the PATs are bounded, allowing the satisfaction of working ranges required from different types of PATs, which has not been considered in available MINLP models. In the case studies, for the first WDS benchmark, we found the same locations of PATs as those reported in the literature, but with a bit higher leakage reduction, i.e., an average of 7.88 L/s compared with 7.72 L/s in the case of three PATs and 6.70 L/s compared with 6.50 L/s in the case of two PATs. For the real WDS in Thainguyen, by installing four PATs in the WDS, we achieve an energy recovery of 1958.0 kWh/day and reduce nearly half of the average excessive pressure. Such results will help water utilities to choose appropriate number of PATs to be installed with expected power generation and leakage reduction. The energy produced reported in this paper represents a reference value for the water utility. To estimate the real energy production, proper types of PATs should be selected, and PAT characteristics and efficiency curves considered. The selection of suitable pumps working as turbines will be studied in our future work. Funding: This research was funded by the science and technology project with grant number B2017-TNA-31.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviation
The following symbols are used in this paper: