Mathematical Program with Vanishing Constraints for Optimal Pressure Control in Water Distribution Systems
Abstract
Optimizing pressure management to reduce water leakage in water distribution systems (WDSs) is one of the major tasks for water utilities. By regulating the operation of pressure reducing valves (PRVs) installed in WDSs, the pressure in WDSs can be kept under control and thus the water leakage amount can be decreased. Mathematically, the problem of pressure management to water leakage reduction can be formulated as a nonlinear optimization program. To make the optimization model proper for practice, the model of PRVs should be accurate and can describe all its operation modes in practice: active, fully opened, and check valve modes. In the literature, the model can be represented either by a nonsmooth equation with low accuracy or by several complicated constraints. This research developed a highly accurate PRV model based on vanishing constraints. The idea comes from the fact that the model equation representing operations of PRVs in active mode will be vanished as PRVs operate in the check valve mode. The formulated mathematical program with vanishing constraints (MPVCs) can be solved efficiently by using the regularization approach. Several WDSs have evaluated the new PRV model which shows that accurate solutions are obtained with less computation time.
1 Introduction
Water loss is inevitable in all water distribution systems (WDSs) worldwide (Lambert 2002; Mosetlhe et al. 2021). In addition, water resources used for drinking water become scarce due to over excessive exploitation and the invasion of saltwater. For this reason, saving drinking water becomes more urgent and, as carried out, will not only preserve water resources for our generation, but also reduce electricity costs for water treatment and water pumping for WDSs. As an example, in India, nonrevenue water (NRW) is relatively high and ranges from 5% to about 70% of total drinking water production (Rajakumar et al. 2020). To decrease the NRW, many approaches can be taken from detecting leaks, localizing pipe bursts/breaks, and rehabilitating WDSs (Puust et al. 2010; Farley and Trow 2003; Lambert 2002). Among these approaches, pressure control is the most efficient for leakage management for existing networks. For new WDSs, other methodologies could be used to reduce pressure in the systems and reduce leakages (Ulanicki et al. 2000; Taha et al. 2020). The reasons are because the water leakage amount at a node, in nature, is proportional to the pressure at the node, and therefore by keeping the pressure under control, the water leakage flow can be decreased. Besides, as compared to the cost of WDS rehabilitation and leakage detection, pressure control can be accomplished with much less cost (Creaco and Walski 2017; Rajakumar et al. 2020). In WDSs, pressure control can be efficiently carried out by operating pressure reducing valves (PRVs) in an optimal fashion (Vairavamoorthy and Lumbers 1998; Dai and Li 2016; Kiliç 2021; and Dai 2021a) where the optimization model is developed and solved to provide pressure set points for PRV controllers. Development of an accurate optimization model for the problem of optimal pressure control is critical. In general, PRVs can operate in one of three discrete modes: active, fully open, and check valve modes. A PRV installed on pipes serves two functions: 1) maintaining the pressure on the downstream side of the PRV at a preset pressure setting regardless of how high the pressure of the upstream nodes is (i.e., active mode); and 2) avoiding reverse flow, as the pressure on the downstream side of the valve exceeds the pressure on its upstream side (i.e., check valve mode) (Simpson 1999). In addition, it is well known that PRVs can be operated in the check valve mode as water demand varies (Sivakumar and Prasad 2015).
Mathematically, optimal pressure management can be transformed into a nonlinear program (NLP) where state variables are flow rates and nodal pressures of WDSs, while decision variables are the pressure settings of PRVs. This optimization problem is also known as the valve setting problem. The constraints of NLPs consist of hydraulic model equations relating to nodal heads and flows of pipes, nodes, and PRVs. Since the NLP model should be accurate and suitable for practice, the formulating model of PRVs should be able to describe two operating functions in active and check valve mode (Dai and Li 2016; Dai 2021b). Many solution approaches, from modeling to solving the NLP, have been proposed in the literature. Genetic algorithm (GA) has been extensively employed for addressing the optimization problems for optimal pressure management (Savić and Walters 1996; Araujo et al. 2006; Awad et al. 2009; Mahdavi et al. 2010; Covelli et al. 2016). In Creaco and Pezzinga (2014), a hybrid optimization algorithm combining GA and a linear program (LP) was used to address the optimal localization and pressure settings of PRVs. In the approach, GA was applied to determine optimal locations of PRVs to be placed, while LP was used to determine optimal pressure settings of PRVs. With this combination approach, the speed of the GA has been shown to be significantly increased. Another approach incorporated priori knowledge of the WDS (i.e., potential links for placing PRVs) in GA for addressing the localization and valve setting problem as discussed in Ali (2014). Besides, a harmony search algorithm was used in De Paola et al. (2017) to address the valve setting problem for minimization of the water leakage flows. The results have shown that better performance of the algorithm as compared with GA was achieved. In Mehdi and Asghar (2019), Particle Swarm Optimization (PSO) was used to determine the optimal pressure settings of PRVs to maximize reliability of the WDS. Differentiable Evolution (DE) was also applied in Cao et al. (2019) to compute optimal pressure settings for real time control of PRVs. To this end, the advantages of heuristic algorithms lie in the fact that they do not need the optimization model to be continuous and differentiable so that the discrete operation modes of PRVs and pumps can be easily included in the optimization model. However, they require an extensive computation time. Recently, Mosetlhe et al. (2021) applied a reinforcement learning for control of pressure in a WDS. At first, a greedy algorithm was applied to identify critical nodes targeted for pressure control, and then a controllerbased reinforcement learning approach was used to determine optimal pressure settings of PRVs. Since the reinforcement learning approach does not use as much data as other machine learning algorithms, the computation time can be decreased, thus enabling realtime control of WDSs.
For mathematical programming, modeling discrete states of control devices such as PRVs is not a trivial task. For this reason, the common equations for PRVs only represent the active operation mode of PRVs. In Vairavamoorthy and Lumbers (1998), a sequential quadratic program (SQP) was applied to address the NLP for optimal pressure management where either excessive pressure or the water leakage amount is minimized. The results have indicated that using the objective function for minimization of excessive pressure in formulation of the NLP results in a better optimal solution. Similarly, the same optimization model is applied to regulate pressure in domestic area meters (DMAs) where both boundary and internal PRVs are used to regulate pressure. From the results, a modulation curve relating flows and pressure settings of each PRV is established and thus allows the realtime control scheme to accomplish for each PRV (Ulanicki et al. 2008). In Wright et al. (2014), the model equation of PRVs in active mode was proposed, where the additional head dissipated by PRVs is introduced and considered as decision variable. Moreover, the method of sequential convex program (SCP) was used to address the NLP, aiming to improve its solution. In the framework of the SCP algorithm, starting from a feasible solution, nonlinear constraints are linearized, and the resulted linear program (LP) is solved to update the current solution. As such, the strategy (linearization, updating solution) is accomplished iteratively until the convergence criteria is satisfied (Wright et al. 2015). Ghaddar et al. (2017) formulated the NLP as a polynomial optimization problem and the approach of semidefinite (SDP) relaxation was applied to solve the optimization problem. Recently, Nerantzis et al. (2020) developed optimization problems for deducing control curves for both variable speed pumps and pressure reducing valves in the WDS. Based on the constructed curves, the control of real networks can be accomplished automatically. For solving the NLP, the linear relaxation with a bound tightening approach was employed to approximate the nonlinear model equations by a set of linear constraints so that a global solution can be found. Numerical experiments have shown that the proposed approach provided the best lower bound of optimal solution as compared with methods used in the BARON solver and water network branchandbound (WNBB). In these works, the check valve mode of PRVs was not considered in the NLP formulation, which may lead to an unsuitable solution. This is because the optimization model will also lead to an unsuitable solution, especially as the check valve mode of PRV occurs. As a result, the solution may be unsuitable and results in high excessive pressure, as compared with the solution of the optimization considering the check valve modes of PRV. Developing a mathematical model of PRVs representing full operating modes of real PRVs is necessary. In this context, Dai and Li (2016) proposed a compact and nonsmooth model equation for PRVs. Recently, the same author (Dai 2021b) proposed a more accurate PRV model based on KarushKuhnTucker conditions for approximate quadratic function. Zhang et al. (2021) introduced a novel DMA sectorization approach to enhance the optimization of background leakage reduction. Additionally, in the work by Jenks et al. (2023), a combination of pressure management and selfcleaning capacity was presented, involving the simultaneous regulation of pressurereducing valves and automatic flushing valves.
In this paper, our objective is to develop a simple and compact model for PRVs. The model is then used for formulation of the optimizing pressure control problem for effective pressure management. It is recognized that as PRVs work in check valve mode (i.e., prevent reverse flow), the model equation for the active mode ought to be vanished. From this recognition, vanishing constraint formulation is applied for modeling operations of PRVs. As a result, the NLP formulated for optimal pressure management is a type of mathematical program with vanishing constraints (MPVCs). This formulation has not been applied in the field of optimal pressure management. Numerical experiments from several case studies have been carried out revealing that the solution of MPVCs is accurate and requires less computation time as compared with the solution of NLPs formulated with the existing PRV models.
The paper is organized as follows: Section 2 presents the formulation of a mathematical program with vanishing constraints (MPVCs) for optimal pressure management, where a new formulation of a PRV model using vanishing constraints is used. A smoothing and regularization approach for solving the MPVCs is given in Section 3. Case studies are outlined in Section 4, and the Conclusion is presented in Section 5.
2 Optimization problem formulation for optimal pressure management
Optimal pressure management is transformed into a nonlinear program (NLP). In general, we consider a WDS with N_{n} nodes, N_{p} links, N_{prv} PRVs, and N_{r} reservoirs. The NLP is formulated to minimize the excessive pressure at all nodes in the time horizon while ensuring that the hydraulic model constraints, as well operational constraints, are satisfied (Germanopoulos and Jowitt 1989). The excessive pressure is defined herein as the deviation of operating pressure and the minimum required to keep at nodes. In addition, as the excessive pressure is reduced, the leakage flow is decreased.
2.1 Objective function
(1) 
In Equation 1, and equal nodal head and its minimum required value, respectively. is defined as the excessive pressure at node i; T is the time horizon (i.e., 24 hours); k is the time interval; and F is the objective function and defined as the total excessive pressures of all demand nodes in the WDS.
2.2 Hydraulic constraints
The continuity equation at node i for time interval k:
(2) 
In Equation 2, at time interval k, d_{i,k }is defined as the demand at node i; Q_{i,j,k} is the flow on pipe i j coming out from node i, while Q_{j,i,k} is the flow on pipe j i coming into node i.
The leakage amount associated with node i, l_{i,k }is a pressure dependent additional demand, and it is mathematically modeled based on the equation of flow through an orifice (Araujo et al. 2006).
(3) 
Where:
i  =  1, ...,N_{n }, 
L_{i,j}  =  length of pipe i,j, 
p_{i}  =  pressure at node i, and 
γ  =  leakage exponent. 
In Equation 3, the coefficient C_{L}, in general, is estimated through the amount of leakage flow at the time of minimum night flow (Araujo et al. 2006).
2.3 Constraints for modeling pipes
The energy equations for pipes connecting node i to node j are:
(4) 
Where the head loss (ΔH_{i,j,k}) through a pipe connecting node i to node j is computed either by the HazenWilliams equation (Bhava and Gupta 2006):
(5) 
or by the DarcyWeisbach equation (Bhave and Gupta 2006):
(6) 
The friction factor f_{i,j} in Equation 7 is implicitly calculated using the ColebrookWhite equation (Bhave and Gupta 2006):
(7) 
Where:
Re  =  Reynolds number, 
D_{i,j}  =  diameter, and 
ε_{i,j }and L_{i,j}  =  roughness and length of the pipe, respectively. 
Equations 5, 6, and 7 are written in metric units.
2.4 Hydraulic model equations for PRVs
Existing models of pressure reducing valves
In the literature, there were some PRV model equations used in the formulation of the NLP for optimal pressure management. Dai and Li (2016) proposed a compact PRV model, which can represent three operation modes (active, open, and check valve) of PRVs in practice. The model is:
(8) 
(9) 
In Equation 8,
R_{i,j}  =  resistance of the PRV, 
K_{i,j}  =  10 = head loss coefficient, 
g  =  gravity constant, 
0 < v_{i,j,k }≤ 1  =  is a factor representing the opening of the PRV at time interval k, and 
D_{i,j}  =  diameter of the PRV. 
Since the PRV model is nonsmooth at a kink, the interiorpoint approximation method proposed in Gopal and Biegler (1999) is applied to smoothen the max term in Equation 8. The smooth model of PRV given in Dai and Li (2016) is shown herein:
(10) 
The shortcomings of Equation 10 are due to the introduction of the approximation parameter τ, which may lead to inaccurate solutions. In addition, the lower bound of decision variable v_{i,j }cannot be set to very small values due to numerical computation.
Another PRV model proposed in Wright et al. (2014) is compact and simpler, but it can only describe two among the three operating modes of PRVs.
(11) 
In Equation 11, instead of introducing a variable representing openings of PRVs (i.e., v_{i,j }), the additional head loss absorbed by PRV (i.e., ) is introduced, and it is an optimization variable. With the same idea, a relaxed PRV model is proposed in Ghaddar et al. (2017) to avoid the introduction of additional variables (i.e., ). In fact, from Equation 11, we have . Since , hence Equation 11 can be relaxed into the following convex inequality one as:
(12) 
Recently, Dai (2021b) proposed an accurate PRV modelbased complementarity constraint. In fact, the nonsmooth term max(0, H_{i,k }– H_{j,k }) is approximated by a set of KKT conditions of the quadratic program. The PRV model is:
(13) 
The last two equality constraints can be treated as complementarity constraints for numerical computations. The variables in Equation 13 are referred to the work in Dai (2021b). The model in Equation 13 requires many additional constraints, thus making the formulated NLP even more complicated. Especially for large scale WDSs, the formulated NLP will have many optimization variables and will be extremely difficult to be solved by standard NLP solvers.
A new model of pressure reducing valves using vanishing constraints
In this paper, our aim is to construct an accurate and compact PRV model equation with a simple model equation. In line with the function of PRVs, as it is placed on links, it performs different functions either in active mode or check valve mode. As seen in Equation 12, the relaxed inequality is used to represent the active mode operations of PRVs.
In check valve mode, the function of the PRV is to prevent the reverse flow through it, so the mathematical PRV model now must force Q_{i,j,k} to be equal to zero, and the inequality of ought to be vanished. To deal with this issue, we applied the formulation of vanishing constraints as follows:
(14) 
The principal operation of a PRV in Equation 15 is that when the PRVs operate in active mode, there are relations, i.e., and . However, in the case that the PRVs work in check valve mode, the constraint is vanished since Q_{i,j,k }= 0. The PRV model constraint in Equation 15 is formulated in such a way that it can be classified into a vanishing constraint (Achtziger et al. 2013).
2.5 Operational constraints
The optimization problem is to minimize excessive pressure while ensuring that the WDS will be operated to the satisfaction of water supply services, i.e., sufficient pressures and flows are guaranteed at the furthest water supply nodes. For this reason, in the NLP, we proposed the following bound constraints on nodal heads and link flows:
Nodal head constraints
(15) 
Pipe flow constraints
(16) 
Pressure reducing valve flow constraints
(17) 
Where H^{L} and H^{U} are the lower and upper bounds of the head at node i, respectively. Similarly, Q^{L} and Q^{U} are the lower and upper bounds of the flows through pipe i j.
2.6 Constant reservoir heads
The head of reservoir i is considered as a constant :
(18) 
Where:
i  =  1, ..., N_{r} 
Since the above NLP consists of vanishing constraints, as shown in Equation 14, the formulated optimization for optimal pressure management is a type of mathematical program with vanishing constraints (MPVCs) (Achtziger et al. 2013). The MPVCs can be written in general as follows:
(19) 
Where h_{i }(x) and g_{i }(x) are constraints for junctions and pipe lines, respectively; G_{i }(x) and H_{i }(x) stands for and Q_{i,j,k}, respectively. In the next section, we will present the regularization approach in Achtziger et al. (2013) to address the MPVCs efficiently.
3 A smoothing regularization solution approach for the MPVCs
Solving the MPVCs in Equation 20 is complicated due to the vanishing constraints. Replacing the vanishing constraint with a smooth function with which the formulated MPVCs can be efficiently solved by standard NLP solvers is necessary. For the sake of simplification, we assign a = G_{i }(x) and b = H_{i }(x). We used a function Ψ:R^{2}→R as proposed by Achtziger et al. (2013) with property:
(20) 
In particular,
(21) 
The function Ψ in Equation 22 is nonsmooth due to the terms of max. and min. In this study, a smoothing and regularization function is developed so that the formulated NLP can be solved by optimization algorithmbased gradient methods. The smoothing form of Ψ ^{t}(a,b) is defined in Equation 22 with a regularized parameter t (Achtziger et al. 2013), and the parameter is chosen as a small positive value, to make Ψ(a,b) smoothed.
(22) 
The function Ψ ^{t }has the following properties (Achtziger et al. 2013):
for all
The gradient is
The following relations are hold
To satisfy the vanishing constraints in Equation 20, the MPVCs are now regularized into a smooth nonlinear program (NLP (t)) by using constraint Ψ ^{t}(a,b) ≤ t as follows:
(23) 
The regularization scheme proposed enlarges the feasible region and the feasible solution set of the program NLP (t) is always nonempty. To this end, to address the MPVCs, we solve a sequence of NLP (t_{k}) where t_{k} is gradually decreased, i.e., t_{k}: = σ_{t }t_{k1, }σ_{t } < 1.0. Inline with this point, the computation framework is given in Table 1. In the table, x^{k} is denoted as the solution of the NLP (t) at the value of t_{k} while ε is chosen as a small positive value.
Table 1 The computational framework for MPVCs.
Data: x^{0} an initial point, t_{0} initial parameters, σ_{t }ε (0,1) parameters update 
1. Set k: = 0, t_{k}: = t_{0} 
2. Repeat 
3. Until t_{k+1} ≤ ε and return x^{k+1} 
4 Case studies
4.1 Case Study 1
In this section, the MPVCs computational framework was applied to solve the problem of optimal pressure management for a real largescale WDS in Thainguyen City, Vietnam, as shown in Figure 1. This WDS was also studied in Dai (2021b) for optimal pressure management consisting of 114 pipes, 6 PRVs, 103 nodes, and 4 reservoirs. The heads of 4 reservoirs at nodes 150, 151, 152, and 153, are assumed to be constant and equal to 111.74 m, 75.21 m, and 69.57 m, respectively. Also, the 24 nominal demand pattern factors are given in Table 2. The formulated MPVCs aim to determine the optimal pressure settings of PRVs to minimize the excessive pressure in the WDS while ensuring that all hydraulic model bound constraints are satisfied. In addition, the minimum allowable pressure at demand nodes is 17.0 m (Dai 2021b). The MPVCs are solved with t_{0 }= 1, σ_{t }= 0.1, and ε = 0.0001. An NLP solver, IPOPT (Wächter and Biegler 2006), is employed to solve the sequence of NLP(t_{k}). All computation experiments are accomplished on computer desktop Intel(R) Core(TM)i56500CPU@3.2GH,RAM 8GB.
Figure 1 A real water distribution system in Thainguyen City.
Table 2 Daily water demand pattern factors.
Time [h]  Demand pattern factors  Time [h]  Demand pattern factors 
1  0.36  13  1.20 
2  0.36  14  1.15 
3  0.36  15  1.15 
4  0.58  16  1.28 
5  0.82  17  1.3 
6  0.86  18  1.34 
7  1.18  19  1.28 
8  1.20  20  1.20 
9  1.25  21  1.12 
10  1.30  22  0.80 
11  1.32  23  0.68 
12  1.34  24  0.46 
Regarding the optimal solution, the optimal pressure settings for PRVs are shown in Figure 2, Figure 3, Figure 4, and Figure 5. At the beginning of the time horizon, the water demand pattern factors are low (i.e., 0.36), and the pressure settings for PRVs, found by MPVCs, are small. For instance, the pressure setting for PRV 1 at the demand factor of 0.36 is around 22.0 m, and for PRV 3, it is near 20.0 m. Similarly, it is 20.0 m for PRV 100, and PRV 6 maintains a pressure of 30.0 m.
Figure 2 Pressure settings of PRV 1.
Figure 3 Pressure settings of PRV 3 and PRV 5.
Figure 4 Pressure settings of PRV 6 and PRV 104.
Figure 5 Pressure settings of PRV 100.
According to the lowpressure settings of PRVs, the average excessive pressure is decreased as shown in Table 3, i.e., it decreases from 32.59 m (i.e., without PRVs) to 4.32 m (i.e., with PRVs). For the PRV working status, at low demand pattern factors, PRVs 5 and 104 are set to work in fully closed mode, while other PRVs are working in active mode. The pressure settings of PRVs gradually increase as the demand pattern grows. This can be seen for PRV 1, PRV 3, PRV 5, PRV 6, and PRV 100, except that PRV 104 is switched to closed operation mode for almost rising demand patterns.
Table 3 Average excessive pressure – Case Study 1.
Time (h)  Average excessive pressure with PRV control (m)  Average excessive pressure without PRV control (m)  Time (h)  Average excessive pressure with PRV control (m)  Average excessive pressure without PRV control (m) 
1  4.36  32.59  13  10.90  20.18 
2  4.36  32.59  14  10.13  21.09 
3  4.36  32.59  15  10.13  21.09 
4  4.49  29.57  16  12.16  18.65 
5  6.35  26.36  17  12.54  18.25 
6  6.73  25.81  18  13.28  17.45 
7  10.59  20.54  19  12.18  18.65 
8  10.90  20.18  20  10.90  20.18 
9  11.73  19.23  21  9.68  21.62 
10  12.54  18.25  22  6.16  26.64 
11  12.90  17.86  23  5.14  28.23 
12  21.57  25.77  24  21.34  28.01 
Since the pressure in the WDS tends to be reduced in high demand patterns, the PRVs are now regulated to not only reduce unnecessarily high pressure, but also guarantee minimum pressure at nodes, especially for the boundary supplying water nodes. For this reason, as in in these Figures, pressure settings found by MPVCs are at reasonable levels, but they are not too low. The maximum pressure settings for PRVs are around 50.0 m.
Besides the control of excessive pressure reduction, the satisfaction of vanishing constraints for PRVs should be verified and, herein, we demonstrate them in Figure 6. It can be seen that vanishing constraints in Equation 15, i.e., G_{i }(x) H_{i }(x) ≤ 0, are satisfied. For the operating status of PRVs, as seen in the figure, PRV 104 operates mostly in fully closed mode from 3:00 am to 23:00 pm. PRV 100 and PRV 5 operate in either closed or active modes in some time intervals. PRV 1, PRV 3, and PRV 6 work in active modes during the entire time horizon.
Figure 6 Vanishing constraints of PRVs.
To demonstrate the efficacy of the MPVC’s formulation, we next compare the optimal solutions obtained from solving the MPVCs (i.e., formulated with our new PRV model) with the ones obtained from solving the NLP formulated with a PRV model as proposed in Dai (2021b). The comparison results are given in Table 4. When using both PRV models, the same optimal solutions are achieved. However, with our PRV modelbased vanishing constraints and with the same setup conditions, IPOPT took much less computational time. This is reasonable because the formulated MPVCs have less constraints than that using the PRV model as seen in Dai (2021b).
Table 4 Comparisons of PRV modeling formulation – Case Study 1.
Regularized parameter  Formulation of MPVCs  PRV model in Dai (2021b)  
t  Objective function value  CPU time [s] 
Objective function value  CPU time [s] 

1  21089.57  2.18  20754.134  3.86  
0.1  21185.14  0.60  20754.136  2.63  
0.01  21173.82  0.64  20754.138  1.03  
0.001  21128.83  1.05  21146.661  3.15  
0.0001  21147.85  0.94  21166.615  2.16 
4.2 Case Study 2
To evaluate the efficacy of the MPVC formulation for optimal pressure management, we use it to determine the optimal pressure settings for a benchmark WDS. The WDS was studied in Vairavamoorthy and Lumbers (1998) and in Araujo et al. (2006) for valve setting problems consisting of 37 links, 25 nodes, and 3 reservoirs, as shown in Figure 7. To control the pressure, 5 PRVs are placed on links: 1, 11, 21, 20, and 29, respectively. The WDS data is tailored in Vairavamoorthy and Lumbers (1998) and Araujo et al. (2006). In addition, the discharge coefficient C_{L} is 10^{5} and the leakage exponent parameter γ is set to 1.18, as used in Araujo et al. (2006). The MPVCs are formulated and solved to determine the optimal pressure settings of PRVs and their corresponding operating states to minimize excessive pressure in the WDS. The excessive pressure is defined as the sum of all deviations between the working pressure (H_{i }) and the minimum allowable pressure (, i.e., 30.0 m) at all nodes. Similarly, IPOPT (Wächter and Biegler 2006) is employed to solve the sequence of NLP(t_{k}). All computation experiments are accomplished on a desktop computer with Intel(R) Core(TM)i56500CPU@3.2GH,RAM 8GB.
Figure 7 A benchmark water distribution system.
The MPVC is formulated with the PRV models using vanishing constraints and the computation is accomplished with t_{0 }= 1, σ_{t = }0.1, and ε = 0.0001. The optimized pressure settings for PRVs are shown in Figures 8 and 9 for PRVs located on links 1, 11, and 21, respectively.
Figure 8 Optimized pressure settings for PRV 1.
Figure 9 Optimized pressure settings for PRV 11 and PRV 21.
During lowdemand patterns, since the pressure tends to be high, the optimized pressure settings for PRVs will be set to low values to decrease more of the excessive pressure in the WDS. In particular, the pressure settings of PRV 1 are nearly 30.0 m for the entire time horizon, while for PRV 11, the pressure settings are set around 30.2 m from 3:00 am to 5:00 am and are set to higher values from 8:00 am to 11:00 am. The higherpressure setting guarantees the minimum pressure satisfaction at boundary nodes. In the same fashion, for PRV 21, the pressure settings are about 32.0 m. Also, according to the optimized result, two PRVs on links 20 and 29 are scheduled to work in fully closed mode for the entire time horizon. The optimized pressure settings for PRVs are regulated according to the variation of water demand patterns. It can be seen in Figure 9 that the pressure settings for PRV 21 are around the value of 32.0 m, while the pressure settings for PRV 11 are changed from 30.0 to 32.0 m. This difference is due to the fact that PRV 11 is near the reservoir and thus its pressure settings are more sensitive to the diurnal change of demand patterns.
It is necessary to evaluate the accuracy of the PRV modelbased vanishing constraints. Herein, we pass the optimized pressure settings, found by the MPVCs, to EPANET 2 (Rossman 2000) as control inputs for carrying out the simulation. The simulated PRV flows will be compared with the optimized ones (i.e., resulting from solving the MPVCs). The comparison results are shown in Figures 10, 11, and 12 revealing that both optimized and simulated flows are almost the same. This indicates that the proposed PRV model based on vanishing constraints is indeed accurate.
Figure 10 Optimized and simulated flows of PRV 1.
Figure 11 Optimized and simulated flows of PRV 11.
Figure 12 Optimized and simulated flows of PRV 21.
Similar to Case Study 1, the satisfaction of vanishing constraints for PRVs are demonstrated in Figure 13. It shows that PRVs on links 20 and 29 are scheduled to operate in fully closed mode, while PRVs on links 1, 11, and 21 are set to active mode with appropriate pressure settings.
Figure 13 Vanishing constraints of PRVs.
For water leakage flow reduction with optimal control operations of PRVs, the hourly water leakage flows are given in Table 5, revealing that the average leakage flow is significantly reduced. In lowdemand patterns, more excessive pressure is eliminated, and as a result, leakage flow is decreased. To this end, by installing 5 PRVs in the WDS, it is possible to save up to 1,821.00 m^{3} of water flow per day.
Table 5 Water leakage flow – Case Study 2.
Time [h]  Leakage flow with PRV control [L/s]  Leakage flow without PRV control [L/s]  Time [h]  Leakage flow with PRV control [L/s]  Leakage flow without PRV control [L/s] 
1  21.35  28.12  13  21.44  26.62 
2  21.35  28.20  14  21.44  26.57 
3  21.30  28.83  15  21.43  26.49 
4  21.31  28.92  16  21.41  26.38 
5  21.31  28.96  17  21.46  25.83 
6  21.32  29.02  18  21.47  25.90 
7  21.46  27.76  19  21.42  26.45 
8  21.46  27.77  20  21.43  26.50 
9  21.68  25.78  21  21.39  26.97 
10  21.65  25.57  22  21.40  27.05 
11  21.58  25.89  23  21.33  27.91 
12  21.57  25.77  24  21.34  28.01 
The accuracy of the new PRV modelbased vanishing constraints is also investigated by comparing the optimal solution resulted by the MPVCs with the one resulted by the MPCCs, where the PRV model is described by a set of constraints from KKT conditions, as proposed in Dai (2021b). The comparison results are given in Table 6. As seen in the table, both the MPVCs in this paper and MPCCs from Dai (2021b) give the same optimal results. However, when using the MPVCs, IPOPT took less computation time. This is suitable since the MPVCs are formulated using less constraints than the MPCCs formulated using the PRV model in Dai (2021b).
Regarding computational time comparisons, we consider both models of PRVs, as in the previous case study. The results of objective functional values and computation times respective to different values of t are given in Table 6. It is seen that, even with a large t value of 1.0, the solution of the MPVCs is close to the optimal one, since the objective function value is nearly the same as those with a smaller value of t. Both PRV models give the same result, but the PRV modelbased on vanishing constraints outperforms in computation time.
Table 6 Comparisons of PRV modeling formulation – Case Study 2.
Regularized parameter  Formulation of MPVCs  PRV model in Dai (2021b)  
t  Objective function value  CPU time [s] 
Objective function value  CPU time [s] 

1  830.91843  0.62  830.9145  1.00  
0.1  830.91932  0.15  830.0992  1.02  
0.01  830.91954  0.24  830.9193  1.05  
0.001  830.91981  0.30  830.9196  1.37  
0.0001  830.91982  0.50  830.9198  1.45 
4.3 Case Study 3
Next, we considered the largest WDS benchmark, which is the EXNET WDS (Farmani et al. 2004) as shown in Figure 14. The system consists of 2465 pipes, 1892 nodes, 2 reservoirs, and 8 PRVs. For pressure regulation, PRVs are placed on links: 3783, 2938, 3046, 5120, 5162, 2699, and 3786. The MPVCs determine the optimal pressure settings for PRVs to minimize the excessive pressure in the WDS. The WDS was also investigated in Eck and Mevissen (2012), and Dai and Li (2016) for optimal PRV localization and settings. The 24 demand pattern factors are given in Table 7. The minimum allowable pressures at nodes are maintained at 8.0 m. In addition, the nonsmooth DarcyWeisbach equation is replaced by a smooth one, as seen in Burgschweiger et al. (2009), thus allowing formulation of continuous MPVCs.
Figure 14 A benchmark water distribution system.
Table 7 Demand pattern factors.
Time [h]  Demand pattern factors  Time [h]  Demand pattern factors 
1  0.851  13  0.886 
2  0.875  14  0.838 
3  0.747  15  0.936 
4  0.743  16  0.941 
5  0.848  17  0.995 
6  0.999  18  0.99 
7  0.744  19  0.93 
8  0.731  20  0.933 
9  0.763  21  1.000 
10  0.856  22  0.971 
11  0.827  23  0.735 
12  0.813  24  0.856 
Similar to Case Study 1, the MPVCs in this study are addressed with t_{0 }= 1, and σ_{t} = 0.1. To solve the sequence of NLP (t_{k}), IPOPT (Wächter and Biegler 2006) is also employed. All computation experiments are also accomplished on a desktop computer with Intel(R) Core(TM)i56500CPU@3.2GH, RAM 8GB.
After solving the MPVCs, the optimal pressure settings of PRVs are given in Figure 15 and Figure 16. It is seen that at the low demand pattern factors, the pressure settings for PRVs are set to low values. This is reasonable because high pressures appear in nature, as the demand consumption is low, and thus the PRVs need to dissipate much more pressure. Reversely, for high demand pattern factors, where the pressure in nature tends to be low, the pressure settings are set to higher values to ensure sufficient pressure at the furthest nodes. It is also seen that PRV 3046 and PRV 5120 operate in active modes for the wholetime horizon, while the rest of the PRVs operate in either active or closed modes.
Figure 15 Optimized pressure settings for PRV 3046 and PRV 5120.
Figure 16 Optimized pressure settings for PRV 5162 and PRV 2699.
Table 8 shows the average excessive pressures when PRVs are installed and optimized, and where no PRV is installed. This shows that we can reduce the average excessive pressure to a maximum of 5.0 m.
Table 8 Average excessive pressure.
Time [h]  Average excessive pressure with PRV control [m]  Average excessive pressure without PRV control [m]  Time [h]  Average excessive pressure with PRV control [m]  Average excessive ressure without PRV control [m]  
1  32.11  37.16  13  30.44  35.21  
2  30.97  35.83  14  32.71  37.85  
3  37.08  42.26  15  27.88  32.23  
4  37.28  42.44  16  27.62  31.93  
5  32.25  37.32  17  24.78  28.43  
6  24.65  28.17  18  24.96  28.77  
7  37.24  42.39  19  28.20  32.60  
8  37.87  42.97  20  28.04  32.42  
9  36.30  41.53  21  24.63  28.10  
10  31.88  36.89  22  25.99  30.02  
11  33.20  38.43  23  37.67  42.79  
12  33.85  39.14  24  31.88  36.89 
Also, to evaluate the accuracy of the new PRV model, we compare optimized PRV flows (i.e., results from solving the MPVCs), with simulated flows where the optimized pressure settings are input to EPANET 2 (Rossman 2000). The comparison results are represented in Figures 17 and 18 for PRVs on links 3046, and 5120, respectively. For other PRVs, similar results are observed. In this case, for the very large WDS, both optimized and simulated flows are almost the same. This demonstrates the accuracy of the PRV model based on vanishing constraints.
Figure 17 Optimized and simulated flows of PRV 3046.
Figure 18 Optimized and simulated flows of PRV 5120.
Like previous case studies, the vanishing constraints are shown in Figure 19 and Figure 20 for PRVs. I can be seen that all vanishing constraints for PRVs are satisfied (i.e., G_{i} (x) H_{i} (x) ≤ 0). For the operation status of PRVs, PRV 3783 as seen in Figure 20, operates mostly in closed mode because the vanishing constraints are equal to zero. In contrast, PRV 2983, PRV 3046, and PRV 5120 all work in active mode due to the negative values of vanishing constraints. The rest of the PRVs work in either closed mode or active mode.
Figure 19 Vanishing constraints of PRVs.
Figure 20 Vanishing constraints of PRVs.
5 Conclusion
Optimal pressure management for water leakage reduction is emergent for water utilities worldwide. This can be accomplished by optimal control of the operation of pressure reducing valves (PRVs) installed in WDSs. In this paper, we have applied novel vanishing constraints for modeling the discrete operation modes of PRVs. As a result, the optimal pressure management problem can be formulated as a mathematical program with vanishing constraints (MPVCs). Numerical experiments have been carried out to determine the optimal pressure settings for PRVs for optimal pressure management for several WDSs. The results have revealed that using the MPVCs formulation, the solution is accurate, and at the same time, the computation time is decreased as compared to approaches using the existing PRV model in the literature. In addition, with the optimized operations of PRVs, the excessive pressure in WDSs is significantly decreased and the leakage flows are reduced. Water resources have become scarce, demanding strict water conservation measures. The effectiveness of optimal pressure management can be further improved by regulating the operations of Pressure Reducing Valves (PRVs) in realtime. Additionally, the energy dissipation resulting from the reduction of excessive pressure should be converted into electrical energy by using pumps working as turbines (PATs). Moreover, the approach of control of both variable speed pumps and PRVs should be addressed to enhance the overall reduction of water leakage.
References
 Achtziger, W., T. Hoheisel, and C. Kanzow. 2013. “A smoothingregularization approach to mathematical programs with vanishing constraints.” Computational Optimization and Applications 55 (3): 733–767.
 Ali, M.E. 2014. “Knowledgebased optimization model for control valve locations in water distribution networks.” Journal of Water Resources Planning and Management 141 (1): 04014,048.
 Araujo, L., H. Ramos, and S. Coelho. 2006. “Pressure control for leakage minimisation in water distribution systems management.” Water Resources Management 20 (1): 133–149.
 Awad, H., Z. Kapelan, and D. Savic. 2009. “Optimal setting of timemodulated pressure reducing valves in water distribution networks using genetic algorithms.” In Integrating Water Systems (CCWI 2009 Conference), 1–3.
 Bhave, P.R., R. Gupta. 2006. Analysis of water distribution networks. Alpha Science International Ltd., Oxford, UK.
 Burgschweiger, J., B. Gnädig, and M.C. Steinbach. 2009. “Optimization models for operative planning in drinking water networks.” Optimization and Engineering 10 (1): 43–73.
 Cao, H., S. Hopfgarten, A. Ostfeld, E. Salomons, and P. Li. 2019. “Simultaneous sensor placement and pressure reducing valve localization for pressure control of water distribution systems.” Water 11 (7): 1352.
 Covelli C., L. Cozzolino, L. Cimorelli, R. Della Morte, and D. Pianese. 2016. “Optimal location and setting of PRVs in WDS for leakage minimization.” Water Resources Management 30 (5): 1803–1817.
 Creaco, E., and G. Pezzinga. 2014. “Multiobjective optimization of pipe replacements and control valve installations for leakage attenuation in water distribution networks.” Journal of Water Resources Planning and Management 141 (3): 04014,059.
 Creaco, E., and T. Walski. 2017. “Economic analysis of pressure control for leakage and pipe burst reduction.” Journal of Water Resources Planning and Management 143 (12): 04017,074.
 Dai, P.D. 2021a. “A new mathematical program with complementarity constraints for optimal localization of pressure reducing valves in water distribution systems.” Applied Water Science 11 (9): 1–16.
 Dai, P.D. 2021b. “Optimal pressure management in water distribution systems using an accurate pressure reducing valve modelbased complementarity constraints.” Water 13 (6): 825.
 Dai, P.D., and P. Li. 2016. ”Optimal pressure regulation in water distribution systems based on an extended model for pressure reducing valves.” Water Resources Management 30 (3): 1239–1254.
 De Paola, F., E. Galdiero, and M. Giugni. 2017. “Location and setting of valves in water distribution networks using a harmony search approach.” Journal of Water Resources Planning and Management 143 (6): 04017,015.
 Eck, B.J., and M. Mevissen. 2012. “Valve placement in water networks: Mixedinteger nonlinear optimization with quadratic pipe friction.” Report No RC25307 (IRE1209014), IBM Research (September).
 Farley, M., and S. Trow. 2003. Losses in water distribution networks. IWA Publishing, London, UK.
 Farmani, S., D.A. Savic, and G.A. Walters. 2004. “"EXNET" benchmark problem for multiobjective optimization of large water systems.” IFAC Workshop, Venice, Italy.
 Germanopoulos, G., and P. Jowitt. 1989. “Leakage reduction by excess pressure minimization in a water supply network”. In ICE Proceedings, Ice Virtual Library 87, 195–214.
 Ghaddar, B., M. Claeys, M. Mevissen, and B.J. Eck. 2017. “Polynomial optimization for water networks: Global solutions for the valve setting problem.” European Journal of Operational Research 261 (2): 450–459.
 Gopal, V., and L.T. Biegler. 1999. “Smoothing methods for complementarity problems in process engineering.” AIChE Journal 45 (7): 15351547.
 Jenks, B., A.J. Ulusoy, F. Pecci, and I. Stoianov. 2023. “Dynamically adaptive networks for integrating optimal pressure management and selfcleaning controls.” Annual Reviews in Control 55, 2023, 486–497.
 Kiliç, R. 2021 “The strategic development for water loss prevention.” Applied Water Science 11 (2): 1–11.
 Lambert, A. 2002. “International report: water losses management and techniques.” Water Science and Technology: Water Supply 2 (4): 1–20.
 Mahdavi, M.M., K. Hosseini, K. Behzadian, A. Ardehsir, and F. Jalilsani. 2010. “Leakage control in water distribution networks by using optimal pressure management: A case study.” In Water Distribution Systems Analysis 2010, 1110–1123.
 Mehdi, D., and A. Asghar. 2019. “Pressure management of largescale water distribution network using optimal location and valve setting.” Water Resources Management 33 (14): 4701–4713.
 Mosetlhe, T.C., Y. Hamam, S. Du, and E. Monacelli. 2021. “Appraising the impact of pressure control on leakage flow in water distribution networks.” Water 13 (19): 2617.
 Nerantzis D., F. Pecci, and I. Stoianov. 2020. “Optimal control of water distribution networks without storage.” European Journal of Operational Research 284 (1): 345354.
 Puust, R., Z. Kapelan, D. Savic, and T. Koppel. 2010. “A review of methods for leakage management in pipe networks.” Urban Water Journal 7 (1): 25–45.
 Rajakumar, A.G., A.A. Cornelio, and M. Mohan Kumar. 2020. “Leak management in district metered areas with internal pressure reducing valves.” Urban Water Journal 17 (8): 714–722.
 Rossman, L. 2000. EPANET 2 User’s Manual. U.S. Environmental Protection Agency, Washington, DC, EPA.600/R00/057.
 Savić, D.A., and G.A. Walters. 1996. “Integration of a model for hydraulic analysis of water distribution networks with an evolution program for pressure regulation.” ComputerAided Civil and Infrastructure Engineering 11 (2): 87–97.
 Simpson, A.R. 1999. “Modeling of pressure regulating devices: The last major problem to be solved in hydraulic simulation.” In WRPMD'99: Preparing for the 21st Century, 1–9.
 Sivakumar, P., and R. Prasad. 2015. “Extended period simulation of pressuredeficient networks using pressure reducing valves.” Water Resources Management 29 (5): 1713–1730.
 Taha, A.W., S. Sharma, R. Lupoja, A.N. Fadhl, M. Haidera, and M. Kennedy. 2020. “Assessment of water losses in distribution networks: Methods, applications, uncertainties, and implications in intermittent supply.” Resources, Conservation and Recycling 152: 104, 515.
 Ulanicki, B., H. AbdelMeguid, P. Bounds, and R. Patel. 2008. “Pressure control in district metering areas with boundary and internal pressure reducing valves.” In 10th International Water Distribution System Analysis Conference, 17–20.
 Ulanicki, B.P., L.M. Bounds, J.P. Rance, and L. Reynolds. 2000. “Open and closed loop pressure control for leakage reduction.” Urban Water 2 (2): 105–114.
 Vairavamoorthy, K., and J. Lumbers. 1998. “Leakage reduction in water distribution systems: optimal valve control.” Journal of Hydraulic Engineering 124 (11): 1146–1154.
 Wächter, A., and L.T. Biegler. 2006. “On the implementation of an interiorpoint filter linesearch algorithm for largescale nonlinear programming.” Mathematical Programming 106 (1): 25–57.
 Wright, R., E. Abraham, P. Parpas, and I. Stoianov. 2015. “Control of water distribution networks with dynamic DMA topology using strictly feasible sequential convex programming.” Water Resources Research 51 (12): 9925–9941.
 Wright, R., I. Stoianov, P. Parpas, K. Henderson, and J. King. 2014. “Adaptive water distribution networks with dynamically reconfigurable topology.” Journal of Hydroinformatics 16 (6): 1280–1301.
 Zhang, T., H. Yao, S. Chu, T. Yu, and Y. Shao. 2021. “Optimized DMA partition to reduce background leakage rate in water distribution networks.” Journal of Water Resources Planning and Management 147 (10), 04021071.