Optimal Pressure Management in Water Distribution Systems Using an Accurate Pressure Reducing Valve Model Based Complementarity Constraints

Abstract

Water loss according to water leakages in water distribution systems (WDSs) is a challenging problem worldwide. An inappropriate operation of the WDS leads to unnecessarily high pressure distribution in the WDS and thus a large amount of water leakage exists. For this reason, optimal pressure management in WDSs through regulating operations of pressure reducing valves (PRVs) is priority for water utilities. The pressure management can be accomplished in a hierarchical control scheme with high level and low level controllers. While the high level controller is responsible for calculating pressure set points for critical nodes, the task of a low level controller is to regulate the pressures at the critical nodes to the set points. The optimal pressure management in the high level controller can be casted into a nonlinear programing problem (NLP) where PRV models are crucial and determine proper operation of the WDS and quality of overall pressure control. PRV models having been used until now either describe two operating modes (active and open modes) or three operating modes (active, open and check valve modes) with parameter dependence. Such models make the formulated NLP unsuitable for the case PRVs work in check valve modes or resulted in inaccurate NLP solution with unexpected operation modes of PRVs, respectively. Therefore, this paper proposes an accurate PRV model based on complementarity constraints. The new PRV model is parameter-less dependence and is capable of describing complete operation modes of PRVs in practice. As a result, the formulated NLP is general and provides accurate NLP solution. The efficiency of our new PRV model is demonstrated on numerous case studies for optimal pressure management of WDSs.

1. Introduction

According to significant population growth, cities are stressed worldwide predominantly. About 1500 million people in global population have increased in the last 20 years, requiring the increase in the demand of cleaning water supply. Moreover, impacts of climate change and drought had led to some water resources to be exhaustedly exploited and thus an increase in the responsibility to adopt more sustainable management of urban water resources [1], where the problem water leakage control systems is vital [1]. Water losses in water distribution systems (WDSs) comprise real losses and apparent losses. The real losses are mainly caused by leakages at network fittings, pipe joints, breaks and/or bursts in pipes while the apparent losses are due to inaccurate meters and/or unauthorized consumption [2,3,4,5]. Over excessive water losses can lead to service interruption, energy losses and may cause contaminant intrusion due to low or even negative pressure conditions [3,4]. It is due to the fact that reducing all water loss components to zero is impossible both in technical and economically viability [5]. Therefore, water loss components should be accurately assessed and prioritized for a reduction [5,6]. To deal with the water loss problem, many strategies have been proposed such as: detection of pipe bursts and locations of seepages, annual rehabilitation and pressure management [2]. Among them, optimal pressure management (control) is one of the most efficient strategies significantly contributing to the decrease of the water leakage amount [6,7]. It is generally accepted that proper management of pressure in WDSs can lead to more reliable operation [8,9,10,11]. The pressure in WDSs is normally adjusted by control devices (i.e., pressure control valves, pressure reducing valves and variable speed pumps) [11,12,13] in kinds of open loop control or closed-loop control schemes. In open-loop control scheme, pressures at outlets of valves (pressure reducing valves—PRVs) are fixed to constant values or changed according to either time schedule (i.e., time modulation method) or their flows (i.e., flow modulation method) [7,13]. On the contrary, in the closed–loop control scheme, pressures at critical nodes are regulated to desired pressure set points by using proportional (P) or proportional–integral (PI) algorithms [14,15,16].

For improving the efficiency of pressure control strategy, pressure control should be accomplished through a hierarchical control scheme with two control layers as seen in Figure 1. The first control layer is considered as a high level controller where optimal pressure management of WDSs can be formulated as a NLP and it is solved by NLP solvers to determine hourly optimal pressure settings for PRVs (or pressures at critical nodes). Since the WDS is complex and is characterized by largely random variations of water demands. Such the optimal pressures are passed to the second control layer (or a low level controller) as pressure set points where the real time control (RTC) scheme with P or PI algorithms can be performed efficiently in a very high sampling rate (i.e., a control time step can be reduced to about 1 s) to regulate the pressures at the critical nodes (or outlet of valves) to the set points [13,14,15,16,17,18,19,20,21]. One of the effective RTC algorithms was discussed in [16,17] where the predicted head deviation (control error) is eliminated instead of the one at the current time.

This paper focuses on formulating and solving the NLP for determination of optimal pressure valve settings (or desired pressure set points at critical) in the high level controller [13,22]. Many solution approaches having been proposed for modeling and solving the NLP. Genetic algorithms (GAs) combined with a hydraulic simulator, EPANET 2 [23], were also used to determine optimal pressure settings of PRVs or on/off working states of valves for optimal pressure regulations in WDSs [24,25,26,27]. Scatter search algorithm, suitable for solving the optimization problem with a limited number of optimization variables, was applied to determine optimal pressure settings for PRVs in [28]. Using this algorithm, a short computation time for solving the optimization problem can be achieved. Creaco and Pezzinga [29] employed a hybrid algorithm combining GA and linear programming (LP) for optimal valve placements and pressure settings where the computation time (i.e., solving the LP instead of the NLP) was significantly reduced. De Paola et al. [30,31] proposed new methodologies, namely harmony search and music-inspired algorithms, to optimize simultaneously localization of PRVs and their pressure settings. Their performances have been demonstrated better than GA. Particle swarm optimization (PSO) [32] and differential evolution (DE) [33] algorithms were also employed to address the optimal valve setting problems, but for small scale WDSs. Although the heuristic solution approaches can solve complex optimization problems, they normally require a huge computation time, especially for large scale WDSs, hence they are appropriate for off-line control.

Besides heuristic algorithms, optimization based gradient methods were also applied to solve the NLP for identifying optimal pressure settings. Such the methods are suitable for near real time control implementation due to less computation time requirement. Early, the method of successive linear programming (LP) proposed in [34,35,36] was applied to solve the NLP where nonlinear constraints representing head loss equations of pipes and valves are iteratively linearized. The LP is then solved by simplex algorithm. The main shortcomings of this approach are due to a burden computation time and low accuracy of the optimal solution. Vairavamoorthy and Lumbers [37] formulated a NLP for optimal pressure regulation to minimize the leakage flow in WDSs, where the sequential quadratic programming (SQP) approach was used to solve the NLP. To enhance the quality of NLP solution, the method of sequential convex programming (SCP) was applied for determination of optimal pressure settings for PRVs for regulating pressures in DMAs [38]. The idea of this method lies in the fact that the nonlinear and non-convex equality constraints are linearized and the resulting convex optimization problem was solved iteratively until convergent condition is reached. Ulanicki et al. [12,13] optimized operations of both boundary and internal PRVs for 24 h to minimize the leakage flows in the DMAs in UK. As a result, the control flow modulation curves for the PRVs relating between the flows and the pressure settings were deduced. Ghaddar et al. [39] formulated the NLP with quadratic constraints (i.e., nonlinear head loss equations of pipes are approximated by quadratic functions). This NLP is then transformed into a polynomial form, which can be efficiently solved by using hierarchy of semidefinite (SDP) relaxations. In the above mentioned optimization approaches, the model of a PRV placed in a pipe has been commonly described by the Hazen–Williams equation [11,12,13,37], which only describes two (i.e., active and open mode) among three operation modes (i.e., active, open and closed (check valve) mode) of PRVs in practice. Therefore this model cannot be used for represent PRV model in the WDS as it operates in the check valve mode. Recently, the formulated NLP in [38] employed an accurate PRV model but with the assumption that PRV only works either in the active or open mode.

It is well-known that an improper operation of PRVs and PCVs may result in instability events and endanger the overall WDS operation [20,21]. For this reason, in the RTC scheme, Galuppini et al. [20] proposed a gain scheduling approach to account and compensate for the nonlinearities (which causes instabilities) in the system while preserving stability and robustness of the closed-loop in a wide operating region. Recently, the same authors presented an efficient control scheme where the PI algorithm is enhanced with a low-filter and a smith predictor to significantly improve the robustness of the control scheme, while, at the same time, help reducing the cost of control. In the high level controller, in order to avoid unexpected operations of PRVs (and thus WDS operation) and low control performance, the mathematical model of PRVs used for formulation of the NLP must be accurate and describe any circumstances of PRV operations occurring in practice. Although Dai and Li. [10] proposed a PRV model, which can describe complete three operation modes of PRVs (i.e., active, open and check valve mode), it is inaccurate. In addition, according to numerical experiments, the lower bounds of decision variables representing the openings of PRVs cannot be set to a very small value (i.e., 1 × 10⁻²⁰), which again contributes to low accuracy of the NLP solution.

In this paper, we propose an accurate PRV model for formulation of the NLP, which overcomes the shortcomings of existing PRV models. In fact, the new PRV model also uses a non-smooth function for describing the three operating modes of PRVs, but instead of using the interior point method [40] for smoothing the non-smooth term, the approximated function is expressed as a solution of a quadratic program (QP). The resulting Karush–Kuhn–Tucker (KKT) conditions of the QP are constraints including complementarity ones, which will be used for formulation of the NLP. Such the approximation avoids the introduction of parameter (as in the interior point method) and therefore it gains the highest accuracy. More importantly, the proposed PRV model is a general model and can be applied to solve any circumstances of PRV configurations in WDSs where the existing PRV model cannot or results in inaccurate solution. The inaccurate solution may lead to unexpected operations of PRVs. Numerical experiments with our new PRV model have demonstrated that the new PRV model outperforms the existing one with better solution.

The remainder of this paper is organized as follows. In Section 2, the existing PRV models are presented. Based on analyzing the shortcomings of these models, the new PRV model is proposed in this section. The formulation of the NLP with the newly proposed PRV model for optimal pressure management is given in Section 3. Three case studies are carried out in Section 4 and conclusions are presented in Section 5.

2. Existing Mathematical Model of PRV and the Newly Proposed PRV Model

2.1. Existing Model of Pressure Reducing Valves

In practice, a PRV has three operating modes, which are active (i.e., the PRV maintains the downstream pressure at the preset value), open (i.e., when the downstream pressure is lower than the pressure setting) and check valve modes (i.e., when upstream pressure is lower than the downstream pressure) [1,10,33]. Consider a PRV placed on link ij, we denote $H_{i,k}$ as the nodal upstream head and $H_{j,k}$ as the nodal downstream head. Several models of PRVs have been used so far for formulation of the NLP [10,12,37,38,39]. One of them is the PRV model proposed in [10] describing three operation modes of PRVs in practice, which is

$$\max \left(0,H_{i,k}-H_{j,k}\right)-\frac{1}{v_{i,j}}{R}_{i,j}{Q}_{i,j,k}^2=0; Q_{i,j,k}\ge 0$$

(1)

$$R_{i,j}=\frac{8K_{i,j}}{{\pi }^{2}g{D}_{i,j}^4}$$

(2)

where $R_{i,j}$ is the resistance of the PRV and $K_{i,j}=10$ is the head loss coefficient of the PRV when it operates at the fully opened mode [10]; $D_{i,j}$ is the diameter of the PRV; g is the gravity constant and $0<v_{i,j}\le 1$ represents the opening of the PRV. To employ the non-smooth Equation (1) for formulation of the NLP, the interior-point approximation method proposed in [40] is applied to smoothen the left hand side of Equation (1). To the end, the model equation representing the PRV is [10]

$$\frac{1}{2}\left(H_i-H_j+{\left({\left(H_i-H_j\right)}^2+{\tau }^2\right)}^{\frac{1}{2}}\right)-\frac{1}{v_{i,j}}{R}_{i,j}{Q}_{i,j}^2=0;Q_{i,j,k}\ge 0$$

(3)

The introduction of a small–positive parameter $\tau$ is necessary in order to guarantee that derivative values of the left hand side of Equation (3) with respective to $H_i$ and $H_j$ exists at a kink (i.e., $H_i=H_j$). Since the model in Equation (3) depends on the value of parameter $\tau$, solving the formulated NLP may result in inaccurate NLP solution and may lead PRVs to be operated in unexpected operation modes. As an example, when $H_i\le H_j$ occurs, flow $Q_{i,j}$ can be non-zero (i.e., PRVs operate in active mode), while in practice, it must be zero (PRVs operate in closed mode). Additionally, according to numerical experiments, the lower bound of $v_{i,j}$ cannot be set to a very small value (i.e., 1 × 10^{− ×}), which once again lowers the accuracy of the NLP solution.

Another PRV model in Equation (4) used in [38] describes two among three operating modes of PRVs.

$$H_{i,k}-H_{j,k}-R_{i,j}{Q}_{i,j,k}^2-{\delta }_{ij}=0 ;Q_{i,j,k}\ge 0$$

(4)

where ${\delta }_{ij}\ge 0$ is the additional head loss absorbed by PRV and is considered as the optimization variable; $R_{i,j}$ is the resistance of the PRV expressed as in Equation (2).

This model equation exactly describes operations of PRV in either active or open modes. Therefore, it cannot be used to represent the PRV model in the formulation of the NLP as PRVs are worked in the check valve mode for preventing reverse flows. However, the advantage of this model is that the upper bound of ${\delta }_{ij}$ can be set to a very large value without causing numerical difficulties for NLP solvers. In our case studies, the upper bound value of ${\delta }_{ij}$ was set to 1 × 10²⁰.

2.2. An Accurate PRV Model Based Complementarity Constraints

Motivated from the shortcomings of existing PRV models, in this section, we propose a new and accurate PRV model using complementarity constraints. In particular, the new PRV model is given in Equation (5)

$$\max \left(0,H_{i,k}-H_{j,k}\right)-R_{i,j}{Q}_{i,j,k}^2-{\delta }_{i,j}=0; Q_{i,j,k}\ge 0;{\delta }_{i,j}\ge 0$$

(5)

Now, our aim was to approximate the non-smooth term $\max \left(0,H_{i,k}-H_{j,k}\right)$ exactly by complementarity constraints. For the sake of simplicity, we omitted index k and we approximated the non-smooth term $\beta =\max \left(0,\Delta H\right)$ with $\Delta H=H_i-H_j$ by smooth ones. The approximation is derived basing on lemma 1 in (7). This lemma has been applied in the optimal control of dynamic hybrid systems [41] and optimal mixing in gas networks [42], but until now it has not been used for developing the PRV model. To approximate the non-smooth term, we considered the following function

$$f\left(\mathbf{z}\right)=\max \left(0,\Delta H\right)-\beta =0;\mathbf{z}={\left(\Delta H,\beta \right)}^T$$

(6)

This function is continuous, convex and piecewise affine (PWA) [42]. The general idea is to express this PWA function as a solution of an appropriate parametric optimization problem. This can be accomplished exactly and leads to mild constraints [41,42].

Lemma 1 ([41], Lemma 9). 

Consider$f_i:{ℝ}^m\to {ℝ}^n$ be a continuous piecewise affine function (PWA) such that every component$f_i:{ℝ}^m\to {ℝ}^n$ is convex for all$i\in \left\{1,\dots ,n\right\}$then

$$\underset{\xi \in {ℝ}^n}{f\left(\mathbf{z}\right)\in \mathrm{argmin}}\left\{\frac{1}{2}{\Vert \xi -\overline{f}\left(\mathbf{z}\right)\Vert }_Q^2: \xi \ge f\left(\mathbf{z}\right)\right\}$$

(7)

for any positive definite n × n diagonal matrix Q and any affine function$\overline{f}:{ℝ}^m\to {ℝ}^n$with$\overline{f}\left(\mathbf{z}\right)\le f\left(\mathbf{z}\right)$for all$\mathbf{z}\in {ℝ}^m$.

The proof of this lemma is presented in [41] (p. 9). Now, we applied lemma 1 to the PWA function in (6) with m = 2 and n = 1. Choosing diagonal value $Q=1$; $\mathbf{z}={\left(z_1,z_2\right)}^T={\left(\Delta H,\beta \right)}^T$ and choose $\overline{f}\left(\mathbf{z}\right)=-z_2-\eta$ for arbitrary $\eta >0$ and $f\left(\mathbf{z}\right)=\max \left(0,\Delta H\right)-\beta$ then we have $\overline{f}\left(\mathbf{z}\right)=-z_2-\eta =-\beta -\eta <-\beta +\max \left(0,\Delta H\right)=f\left(\mathbf{z}\right)$.

To solve the quadratic program (QP), by ignoring constant term $\frac{1}{2}{\left(\beta +\eta \right)}^2$ in the cost function in (7), the QP becomes

$$f\left(\Delta H,\beta \right)\in \underset{\xi \in ℝ}{\mathrm{argmin}}\left\{\frac{1}{2}{\xi }^2+\left(\beta +\eta \right)\xi : \xi \ge -\beta ,\xi \ge \Delta H-\beta \right\}$$

(8)

The unique solution of the QP is characterized by its Karush–Kuhn–Tucker (KKT) conditions

$$\begin{array}{r} \hspace{1em}\xi +\beta +\eta -{\lambda }_1-{\lambda }_2=0,\\ \beta +\xi \ge 0,\hspace{1em}\beta +\xi -\Delta H\ge 0,\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} {\lambda }_1,{\lambda }_2\ge 0,\\ {\lambda }_1\left(\beta +\xi \right)=0,\hspace{1em}{\lambda }_2\left(\beta +\xi -\Delta H\right)=0.\end{array}$$

(9)

Using that $f\left(\mathbf{z}\right)=0$ provides $\xi =0$, we obtained another complementarity constraint

$$\begin{array}{r} \hspace{1em}\beta +\eta -{\lambda }_1-{\lambda }_2=0,\\ \beta \ge 0,\hspace{1em}\beta -\Delta H\ge 0,\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} {\lambda }_1,{\lambda }_2\ge 0,\\ {\lambda }_1\beta =0,\hspace{1em}{\lambda }_2\left(\beta -\Delta H\right)=0.\end{array}$$

(10)

where ${\lambda }_1$ and ${\lambda }_2$ are internal multipliers. It is well-known that solving such the equality constraints is difficult, therefore we equivalently replaced them by the following inequality ones

$${\lambda }_1\beta \le 0,\hspace{1em}{\lambda }_2\left(\beta -\Delta H\right)\le 0$$

(11)

To the end, the new PRV model based complementarity constraints is

$$\begin{array}{r}\beta -R_{i,j}{Q}_{i,j}^2-{\delta }_{i,j}=0,\\ \beta +\eta -{\lambda }_1-{\lambda }_2=0,\\ \beta \ge 0,\hspace{1em}\beta -\Delta H\ge 0,\\ {\lambda }_1,{\lambda }_2\ge 0;\eta >0;{\delta }_{i,j}\ge 0,\\ {\lambda }_1\beta \le 0,\hspace{1em}{\lambda }_2\left(\beta -\Delta H\right)\le 0.\end{array}$$

(12)

As compared with PRV model in (3), our proposed PRV model equation is parameter less and polynomial, which is beneficial for NLP solvers. In addition, we use a decision variable ${\delta }_{i,j}$ like in [38] (instead of $v_{i,j}$) for representing additional head absorbed by PRVs that the NLP solver can handle easily. It can be concluded that with our new PRV model, we expect that the optimal solution will be accurate and unexpected operations of PRVs can be avoided. For numerical experiments, the lower bound value of arbitrary variable ${\eta }_{i,j,k}$ is chosen as 1 × 10⁻³.

3. Problem Formulation for Optimal Pressure Management

The aim of this section was to formulate the NLP with our newly proposed PRV model where the total excessive pressure was minimized while model equality and inequality constraints describing the WDS model and operational constraints were satisfied. Besides, we also formulated the optimization problem with existing PRV models for purposes of comparisons.

3.1. Objective Function

The objective function to be minimized is defined as the excessive pressure at all nodes in the WDS in the optimization time horizon T [10,35].

$$\underset{\mathbf{H},\mathbf{Q}}{\min }F={\displaystyle \sum _{i=1}^{NJ}{\displaystyle \sum _{k=1}^T\left(H_{i,k}-H_{i,k}^L\right)}}$$

(13)

where $NJ$ is the total number of nodes and T = 24 h is the time horizon. We considered a WDS with $NP$ pipes, $NJ$ junction nodes, $NR$ reservoirs and $NPRV$ pressure reducing valves (PRVs); $H_{i,k}$ is head at node i, at time interval k and $H_{i,k}^L$ is the minimum allowable head at node i.

3.2. Constraints

3.2.1. The Continuity Equation at Node i

$${\displaystyle \sum _{j,k}{Q}_{j,i,k}-d_{i,k}-l_{i,k}=0};i=1,\dots ,NJ$$

(14)

The leakage amount, $l_{i,k}$, associated to node i [8].

$$l_{i,k}=C_L{L}_{t,i}{p}_{i,k}^{\gamma }$$

(15)

$$p_{i,k}=H_{i,k}-E_i;L_{t,i}=0.5{\displaystyle \sum _j{L}_{i,j}}$$

(16)

where $p_{i,k}$ and $E_i$ are the pressure and elevation of node i, respectively; $L_{i,j}$ is the length of the pipe connecting node i to node j; $Q_{j,i,k}$ is the flow rate from node j to node i at time interval k; $d_{i,k}$ is the demand at node i, at time interval k; $L_{t,i}$ is the total length of links connected to node i; $p_{i,k}=H_{i,k}-E_i$ and $H_{i,k}$ is the static pressure and total head at node i, respectively and $\gamma =1.18$ is the leakage exponent, $C_L$ is the discharge coefficient of the orifice and has the unit of $l/\left(s.m^{1+\gamma }\right)$.

3.2.2. The Energy Equation for the Pipe Connecting Node i to Node j

$$H_{i,k}-H_{j,k}-\Delta H_{i,j,k}=0; i,j=1,\dots ,NP$$

(17)

where $H_{i,k}$ is the head at node i; $\Delta H_{i,j}$ is the head loss, which can be computed by the Hazen–Williams equation

$$\Delta H_{i,j.k}=\frac{10.67L_{i,j}}{D_{i,j}^{4.87}{C}_{i,j}^{1.852}}{Q}_{i,j,k}{\left|Q_{i,j,k}\right|}^{0.852}$$

(18)

In Equations (6) and (7) $L_{i,j}$, $D_{i,j}$ and $C_{i,j}$ are the length, diameter and the Hazen–Williams coefficient of the pipe ij, respectively.

3.2.3. Model Constraints for PRVs

In Section 2, the model of PRV is represented by

$$\begin{array}{l}{\beta }_{i,j,k}-R_{i,j}{Q}_{i,j,k}^2-{\delta }_{i,j,k}=0\\ {\beta }_{i,j,k}+{\eta }_{i,j,k}-{\lambda }_{1i,j,k}-{\lambda }_{2i,j,k}=0\\ {\beta }_{i,j,k}\ge 0,\hspace{1em}{\beta }_{i,j,k}-\Delta H_{i,j,k}\ge 0\\ \Delta H_{i,j,k}=H_{i,k}-H_{j,k}\\ {\lambda }_{1i,j}{\beta }_{i,j}\le \rho ,\hspace{1em}{\lambda }_{2i,j}\left({\beta }_{i,j}-\Delta H_{i,j}\right)\le \rho \\ i,j=1,\dots ,NPRV\end{array}$$

(19)

3.2.4. Bound Constraints for Flows and Heads

$$H_i^L\le H_{i,k}\le H_i^U$$

(20)

$$Q_{i,j}^L\le Q_{i,j,k}\le Q_{i,j}^U$$

(21)

$$0\le Q_{i,j,k}\le Q_{i,j}^U; i,j=1,\dots ,NPRV$$

(22)

$$0\le {\delta }_{i,j,k},{\lambda }_{1i,j,k},{\lambda }_{2i,j,k}; i,j=1,\dots ,NPRV$$

(23)

$$0<{\eta }_{i,j,k} ;i,j=1,\dots ,NPRV$$

(24)

3.2.5. The Reservoir Water Levels

$$H_{i,k}=\overline{H}; i=1,\dots ,NR$$

(25)

In the above NLP, optimization variables are flow rates ($Q_{i,j,k}$), nodal heads ($H_{i,k}$), additional head loss absorbed by PRVs (${\delta }_{i,j,k}$) and internal multipliers (${\delta }_{i,j,k},{\lambda }_{1i,j,k},{\lambda }_{2i,j,k},{\eta }_{i,j,k}$). Among them, the decision variable is ${\delta }_{i,j,k}$ and the rest are state variables.

To efficiently solve such an NLP, the computational framework for solving the mathematical program with complementarity constraints (MPCCs) was applied; in fact, complementarity constraints ${\lambda }_{1i,j,k}{\beta }_{i,j,k}\le 0,\hspace{1em}{\lambda }_{2i,j,k}\left({\beta }_{i,j,k}-\Delta H_{i,j,k}\right)\le 0$ were regularized to ${\lambda }_{1i,j,k}{\beta }_{i,j,k}\le \rho ,\hspace{1em}{\lambda }_{2i,j,k}\left({\beta }_{i,j,k}-\Delta H_{i,j,k}\right)\le \rho$ with gradually decreased values of $\rho$; in particular $\rho =1.0$; $\rho =0.01$; $\rho =0.001$; $\rho =0$ was used. As a result, a sequence of the NLPs regularized with these values of $\rho$ was solved. The initial guest of the current NLP was taken from the solution of the previous one. This solution approach ensures that a better solution is found and that the exact PRV model equation is satisfied at a stationary point (i.e., solution of the NLP with $\rho =0$).

4. Case Studies

4.1. Case Study 1: Optimal Pressure Management for an Illustrative Water Distribution System with Multi Reservoirs

In this case study, we considered an illustrative WDS consisting of 14 pipes, 4 PRVs, 16 nodes and 3 reservoirs as seen in Figure 2. The purpose of considering this simple WDS is to evaluate the performance of our proposed PRV model, which can completely describe behaviors of PRVs in practice while providing a fair NLP solution. All pipes had the Hazen–William coefficient of 100.0 and all junction nodes had their elevations of 20.0 m. This WDS was simple, but PRVs were connected in a complicated configuration. The data of links and nodes are given in

Table 1. Data of pipes and nodes of the water distribution system (WDS).

Pipe ID	Start Node	End Node	Length (m)	Diameter (mm)	Node ID	Demand (L/s)
1	1	2	100	300	2	100
2	13	12	100	350	3	0
3	4	5	150	200	5	0
4	5	6	500	200	6	300
5	6	3	0.01	200	7	300
6	3	7	200	300	8	100
7	2	3	150	300	9	50
8	12	9	100	300	10	50
9	12	10	200	300	11	30
10	10	11	100	300	12	50
11	11	7	50	300	14	20
12	6	8	300	400	15	10
13	9	8	150	300	1	Source node with total head of 120.00 m
14	14	5	50	300	4	Source node with total head of 100.00 m
15	15	9	50	300	13	Source node with total head of 120.00 m
16	14	15	500	300
17	16	7	250	300

. The 24 demand patterns were assumed as the same as in [8] and are given in

Table 2. Daily demand patterns.

Time (Hours)	Demand Pattern	Time (Hours)	Demand Pattern
1	0.61	13	0.92
2	0.61	14	0.92
3	0.41	15	0.92
4	0.41	16	0.92
5	0.41	17	1.03
6	0.41	18	1.03
7	0.81	19	0.92
8	0.81	20	0.920
9	1.23	21	0.82
10	1.23	22	0.82
11	1.13	23	0.61
12	1.13	24	0.61

. The NLP was formulated to minimize excessive pressure in the WDS while ensuring that the required minimum pressures at all demand nodes were 30.0 m. The resulting NLP model had 1296 continuous variables and 1200 nonlinear constraints.

IPOPT (https://coin-or.github.io/Ipopt/index.html, accessed on 24 February 2021), a popular nonlinear program solver [43], was employed to solve the NLP for optimal pressure management. Additionally, CONOPT4 solver (http://www.conopt.com/, accessed on 24 February 2021) [44] in the General algebraic modeling system (GAMS) [45] was also used to solve the NLP to ensure a good solution is found. All computation experiments are accomplished on computer desktop CPU- Intel (R) Core (TM) i5 3.2 GH, 8 GB RAM. To compare the performance of our new PRV model with the existing ones, we also solved the NLPs formulated with existing PRV models in [10,38].

For solving the optimization problem formulated with our new PRV model, IPOPT took about 12.0 s and the resulting objective function value was 8693.65 m. For optimal operations of PRVs, PRV 7 was set to be closed in the whole time horizon while PRV 4, 13 and 17 were operated in either active or fully closed modes. The pressure settings were given in

Table 3. Optimal pressure settings of pressure reducing valves (PRVs).

Time (Hours)	PRV 4	PRV 13	PRV 17	Time (Hours)	PRV 4	PRV 13	PRV 17
1	Closed	30.00	Closed	13	30.01	30.00	Closed
2	Closed	30.00	Closed	14	30.01	30.00	Closed
3	Closed	Closed	Closed	15	30.01	30.00	Closed
4	Closed	Closed	Closed	16	30.01	30.00	Closed
5	Closed	Closed	Closed	17	30.00	30.05	Closed
6	Closed	Closed	Closed	18	30.00	30.05	Closed
7	30.00	30.30	Closed	19	30.01	30.00	Closed
8	30.00	30.30	Closed	20	30.01	30.00	Closed
9	30.00	31.30	30.18	21	30.00	30.28	Closed
10	30.00	31.30	30.18	22	30.00	30.28	Closed
11	30.00	30.46	30.19	23	Closed	30.00	Closed
12	30.00	30.46	30.19	24	Closed	30.00	Closed

for these PRVs, respectively.

It is seen in the table that, at low demand pattern values, the pressure settings were set to the minimum value of 30.00 m, while they were a bit higher at larger demand pattern values (i.e., time intervals 9 and 10 with PRV 13 and PRV 17). Moreover, to evaluate the accuracy of our proposed PRV model, the optimal pressure settings of PRVs were passed to EPANET 2 (https://www.epa.gov/water-research/epanet, accessed on 24 February 2021) [23] for carrying out simulation of the WDS. The comparisons of PRV flows, i.e., between ones obtained from solving the NLP and ones from simulation of the WDS using EPANET 2, were accomplished and shown from Figure 3a to Figure 3c. It can be seen in these figures that by using our PRV model very high accuracy was achieved. Moreover, the very small discrepancy between the optimized objective function value (8693.65 m) and the simulated one (8694.83 m) indicating that very high accuracy of nodal heads was also attained. At optimal solution, PRV 17 was operated in check valve modes from 1:00 to 8:00 and from 13:00 to 24:00 (i.e., $H_i<H_j$).

Now, we employed the existing PRV model [38] in the formulation of the NLP. According to numerical experiments, we found that in this case, IPOPT failed to solve the NLP, while CONOPT4 in GAMS could solve it but converged to a very bad local solution with the objective function value of 10118.95 m as seen in

Table 4. Comparison of optimal solutions.

	PRV Model in [10]	PRV Model in [38]	PRV Model Based Complementarity Constraints
$\mathbf{Lower} \mathbf{Bound} \mathbf{of} {\mathit{v}}_{\mathit{i},\mathit{j}}$	Objective Function Values (m)	Objective Function Value (m)	Objective Function Value (m)
1 × 10−6	8731.49	10,118.95 *	8694.83
1 × 10−7	8704.98

* by CONOPT4.

.

For the PRV model in [10], although IPOPT can solve the formulated NLP, it results in lower accuracy of the solution. In particular, the resulting objective function value was 8705.48 m larger than the one found by using our new PRV model, which was 8693.65 m. In addition, if we set the lower bound of $v_{i,j}$ to a very small value, namely less than 1 × 10⁻⁸, due to numerical implementation, IPOPT failed to solve the formulated NLP. The same result is observed by CONOPT 4. More importantly, for the case the optimal solution was found (i.e., lower bound of $v_{i,j}$ was larger than 1 × 10⁻⁷), due to the inaccurate solution, PRV 7 was set to work in the active mode, while, in fact, it was set to be fully closed due to optimal solution resulted by using our new PRV model. The comparisons of optimal solutions with different lower bound values of $v_{i,j}$ are given in Table 4. The inaccurate flows through PRV 7 are shown in Figure 3d. The comparisons of cumulative pressure distributions are given in Figure 4. It is seen that the solution of the NLP formulated with the PRV model in [38] had the worst quality.

4.2. Optimal Pressure Management for a Benchmark Water Distribution System

To verify the accuracy and efficiency of our new PRV model, we considered a benchmark WDS comprising of 32 links, 5 PRVs, 22 nodes and 3 reservoirs as shown in Figure 5. This WDS was also studied in [8,10,34,35] for optimal valve setting problems and studied in [45] for optimal pumps as turbines localization. Note that directions for placing PRVs were from node 23 to node 1 for PRV 1; from node 13 to node 12 for PRV 11; from node 12 to node 15 for PRV 21; from node 22 to node 15 for PRV 20 and from node 22 to node 21 for PRV 29.

The data of links and nodes are given in [8]. In addition, the discharge coefficient $C_L$ was ${10}^{-5}$ and the leakage exponent parameter $\gamma$ was 1.18 as used in [8]. The minimum allowable pressures at nodes were maintained at 30.00 m. In this paper, optimal pressure settings of PRVs were determined for 24 demand patterns and reservoir heads. The reservoir heads for 24 h were taken in [8,10] while the 24 demand patterns are given in Table 1. The resulting NLP model had 2280 variables and 2160 constraints. Similarly, IPOPT was used to solve such the sequence of NLPs and it took about 1.82 s. The results of optimal pressure settings of PRV 1, PRV 11 and PRV 21 are demonstrated in Figure 6a while the PRV 20 and PRV 29 were fully closed in the whole time horizon.

Similarly, in order to evaluate the accuracy of our proposed PRV model, the optimal pressure settings of PRVs were passed to the EPANET 2 [23] as control input for carrying out simulation of the WDS again. The results of simulated PRV flows were compared with the optimized ones (i.e., obtained by solving the NLP) and they are shown in Figure 6b and Figure 7a,b for three PRV 1, PRV 11 and PRV 21, respectively. It can be seen in these figures that our new PRV model was indeed accurate. The high accuracy of the NLP solution comes from the fact that our PRV model did not require a parameter $\tau$ as in [10], and in addition, the problem of a very small lower bound value of $v_{i,j}$ was avoided.

With our new PRV model, the resulting objective function value was 830.91 m, while the simulated one (i.e., obtained by simulation using EPANET 2) was 830.99. This demonstrated fair accuracy of nodal heads. It is also seen in Figure 6a that, at low values of demand patterns, the optimized pressure settings of PRVs were small to restrict the pressures at downstream nodes and thus regulated the system pressure. This is reasonable because at the low demand patterns, the high pressure will be appeared in the system if PRVs are not installed. On the contrary, at high values of demand patterns, the system pressure, i.e., in case of no PRVs, tended to be decreased as compared with the case of the low demand pattern values, therefore the optimized pressure settings of PRVs obtained higher values to ensure that minimum pressure at all demand nodes were at least of 30.0 m.

The water leakage flows for 24 h are shown in Figure 8a. It can be seen that much more water leakage flow was decreased as pressure settings of PRVs were optimized. More specifically, with 5 PRVs to be installed, the total water leakage flow reduced from 56,282.00 to 44,461.00 m³/day. This means that we saved an amount of 11,821.00 m³/day as PRVs were installed and operated at appropriate pressures. The decrease of excessive pressure is represented in the cumulative pressure distribution as demonstrated in Figure 8b, as seen, with 5 PRVs, nearly all demand nodes had their pressures lower than 35.0 m, while only about 10 percent of nodes had such pressure when no pressure control (i.e., no PRV was installed) was made.

Now, the PRV model in [10,38] were used for formulation of the NLP. IPOPT was also used to solve these NLPs. The comparisons of optimal solutions are given in

Table 5. Comparisons of optimal solutions.

	PRV Model in [10]	PRV Model in [38]	PRV Model Based Complementarity Constraints
$\mathbf{Lower} \mathbf{Bound} \mathbf{of} {\mathit{v}}_{\mathit{i},\mathit{j}}$	Objective Function Values (m)	Objective Function Value (m)	Objective Function Value (m)
1 × 10−6	834.10	830.99	830.99
1 × 10−7	832.66

. It is seen in the table that our new PRV model results the same solution as the PRV model in [38] while it slightly outperformed the one in [10]. Similarly to case study 1, CONOPT4 in GAMs also failed to solve the NLP formulated with the lower bound values of $v_{i,j}$ less than 1 × 10⁻⁸.

4.3. Optimal Pressure Management for a Large Scale Water Distribution System in Vietnam

In order to demonstrate the efficiency of our new PRV model, we applied it for optimal pressure management of a large scale and real-world WDS in Thainguyen City in Vietnam as shown in Figure 9. This WDS, considered in [46] for study of optimal pumps as turbines localization, consists of 110 pipes, 6 PRVs, 99 nodes and 4 reservoirs. The constant heads of four reservoirs 150, 151, 152 and 153 were 111.74 m, 75.21 m, 69.57 m and 69.46 m, respectively while the 24 demand pattern factors are given in

Table 6. Daily demand pattern factors.

Time (Hours)	Demand Pattern Factors	Time (Hours)	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.30
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

[46]. The EPANET file of this WDS is given in the supplementary part. The optimization problem aims to determine optimal pressure settings of PRVs so as to decrease excessive pressure in the system while ensuring that the minimum allowable pressures at demand nodes are maintained at 17.00 m. The resulting NLP model had 6072 continuous variables and 5928 nonlinear constraints, respectively.

IPOPT [43] was used to solve the NLP and it took about 12.3 s to solve the sequence of NLPs. The resulting objective function value was 21,071.18 m, while the simulated one by EPANET 2 was 20,656.02 m. Such a small difference revealed that a high accuracy of nodal head was attained. The results of optimal pressure settings are shown in Figure 10a,b for 6 PRVs, respectively.

It is observed in these figures that the optimal pressure settings of PRVs changed according to the daily changes of water demand patterns. In general, at low demand pattern values, the pressure settings of PRVs were usually small while they were higher at high demand pattern values. Similarly, to verify the accuracy of the proposed PRV model, we compared the optimized PRV flows (i.e., obtained from the solving the NLP) with the ones obtained by simulation of the WDS using EPANET 2 where the optimized pressure settings of PRVs were used as control input. The comparisons are shown in Figure 11, Figure 12 and Figure 13 for 6 PRVs, respectively. Once again, it is seen that both optimized flows and simulated flows through PRVs were almost the same demonstrating high accuracy of the proposed PRV model.

To demonstrate the decrease of excessive pressure, the cumulative distributions of pressures of all nodes over 24 h for three PRV models are given in Figure 14, respectively. It is seen that by setting pressures of PRVs to appropriate values, a significant excessive pressure amount was reduced. In particular, with 6 PRVs installed and their operations were optimized, nearly 80% of nodes had pressures lower than 30.0 m comparing with about only 10% of nodes had such a pressure value when no PRV was installed. This means that optimal pressure management through regulating operations of PRVs was indeed efficient. The nice performance of our proposed PRV model was also demonstrated in the figure where the cumulative pressure distribution with the range of dominated pressures from 20.00 to 40.00 m had a slightly larger percent than those resulting from other PRV models.

For comparisons with existing PRV models, we employed the PRV model in [10] and the one in [38] for optimal pressure management; the results of both optimal solutions were compared with the one found by our new PRV model and given in

Table 7. Comparisons of optimal solutions.

	PRV Model in [10]	PRV Model in [38]	PRV Model Based Complementarity Constraints
$\mathbf{Lower} \mathbf{Bound} \mathbf{of} {\mathit{v}}_{\mathit{i},\mathit{j}}$	Objective Function Values (m)	Objective Function Value (m)	Objective Function Value (m)
1 × 10−6	22,007.48	20,828.76	20,656.02
1 × 10−7	21,108.71

. It can be seen in the table that using our new PRV model, a better solution with a lower objective function value was achieved.

5. Conclusions

Pressure control can be efficiently accomplished in a hierarchical control scheme with high and low level controllers. In a high level controller, this paper proposed a new mathematical model for PRVs, which is not only capable of describing complete operation modes of PRVs but also provides accurate solution for optimal pressure managing problem. In fact, the non-smooth term in the new PRV model is expressed as a unique solution of a quadratic program. The KKT conditions of the quadratic program, which are constraints consisting of complementarity ones, will be used for formulation of the PRV model and for formulation of the NLP. The important aspect of the new PRV model lies in the fact that its model constraints were in polynomial kinds, which are beneficial for NLP solvers. In addition, it provides an accurate solution and thus avoids unexpected operations of PRVs. Numerical experiments were verified on three case studies revealing that better solutions could be achieved with our new PRV model. In particular, for an illustrative WDS, although PRVs were connected in a complicated configuration, solving the NLP formulated with our new PRV still gave a better solution as compared with the solution of the existing ones. For a benchmark WDS, by optimizing operations of PRVs, it is possible to save an amount of 11821.0 m³ water per day. Finally, for a large scale and real-world WDS, once again, solving the NLP with our new PRV model results in a better solution and a significant excessive pressure reduction was reached.