The effect of downstream valve closure scheduling was analyzed to find optimal closure parameters that lead to the minimal maximum pressure head in the water distribution system. Several valve closure strategies were explored, combining the known valve performance curve (change in flow as a function of change in valve's opened area) with unknown valve closure curve (change in valve's opened area as a function of time). Second-order polynomial curve, power function curve, and piece-wise linear curve were implemented and compared. Genetic algorithm and quasi-Newton (QN) optimization methods were applied. The methodology was tested for three networks, including looped gravitational and pressurized networks. The results demonstrate that flexible multi-parametric valve closure curve and QN optimization method are more effective in minimizing the maximum pressure head in the system.

INTRODUCTION

A water distribution system (WDS) conveys water from the water sources to water consumers. The operation of the WDS is maintained by controlling network active elements (actuators), i.e., pumps and valves, while keeping desired performance levels at network passive elements, e.g., pipes and tanks. Valves can be used to control flow and pressure (e.g., to protect pipes from excessive pressures) or completely block flow (prevent back-flow). Pumps are used to add hydraulic energy to the system (Merritt 1967). The operation of network actuators induces transients in the system in the form of pressure waves traveling from the disturbance thorough the system (Boulos et al. 2005). Improper operation of actuators, e.g., rapid valve closure or pump start-up, can result in severe transients causing significant pressure oscillations (Azoury et al. 1986).

Valve-closure induced transients are among the most destructive transient accidents, for example, Oigawa hydropower station accident in 1950 in Japan (Bergant et al. 2006), San Onofre Nuclear Generation Site Unit 1 accident in 1985 (Tian et al. 2008), a water hammer event at Hanford, in a Department of Energy facility, in 1993 (ORPS Reports RL-WHC-WH-C300EM-1993-0022). Valve manufacturers’ manuals and guidelines advise slow and gradual closure, and lack an analytical or system-defined approach (e.g., Operation and Maintenance Manual for Gate, Globe and Check Valves). Potential damage risks of the propagating wave drove transient researchers to study this issue.

Streeter (1966) was the first to study flow adjustment in a controlled manner during valve closure to ensure that the pressure head in the system remains within predetermined limits. Azoury et al. (1986) investigated the effect of valve-closure schedule on water hammer under turbulent friction conditions for a single pipeline discharging into free air. The authors provided a chart to determine the valve schedule that should be used to yield minimum water hammer pressure. An effective technique for valve-motion synthesis was developed and applied by Goldberg & Karr (1987). The method was tested on several single pipe example problems. Tian et al. (2008) showed that severe slamming between the valve disk and valve seat may occur during the alternate startup of parallel pumps. The authors suggest adapting damping torque to slow down valve closing speed and propose optimum system design to mitigate the potential damage caused by valve-induced water hammer. An optimal closing rule curve was also developed by Bazargan-Lari et al. (2013). The authors used a multi-objective optimization model and Bayesian networks for controlling water pressure during and after valve closure. In this study, non-dominated sorting genetic algorithm-II was used to develop a Pareto front among three objectives: maximum and minimum water pressures and the amount of water passing through the valve during the valve-closing process. The method of characteristics was used to simulate transient flow in all of the mentioned works.

In this study, we assume valve closure time is constant and short, i.e., we consider the case where valve closure can pose potential danger. In practice, this refers to the cases where valve closure time cannot be increased to prevent water hammer. At fixed closing time, valve closing curve choice plays a substantial role in controlling pressure variations in a pipeline (Bazargan-Lari et al. 2013). The goal and contribution of this study is to investigate different closure curves for valve operation and different optimization algorithms for their applicability to transient calculations. Three types of curves were implemented to model valve closure: (1) quadratic function, (2) power function, and (3) piece-wise linear. The optimal operation of valves or pumps was formulated as an optimization problem with the maximum pressure head as the objective function, transient equations as constraints, and operation parameters as decision variables. Two solution approaches were investigated to find optimal parameters for these curves: (1) genetic algorithm (GA) and (2) quasi-Newton (QN) approaches. The following sections describe the suggested approaches and demonstrate their performance using three water networks.

METHODOLOGY

The approach involves four main stages: (1) selection and definition of valve characteristic curve, closure duration and type; (2) hydraulic transient modeling and systems response to valve operation; (3) formulation and initialization of the optimization algorithms, i.e., parameterization of GA and QN; and (4) interpreting the solution results. The approach overview is presented in Figure 1 and highlights the main phases of the approach and the connections between them.
Figure 1

Methodology overview.

Figure 1

Methodology overview.

Transient modeling

A conceptual description of the water hammer phenomenon and the derivation of the mathematical model as a hyperbolic partial differential equation has been developed by Streeter & Wylie (1983), Chaudhry (1979), and Larock et al. (2000) and is commonly applied in research and practice. The one-dimensional unsteady pressure flow equations describing the water hammer effect can be formulated as: 
formula
1
 
formula
2
where V is the average velocity at a cross section, t is time, ρ is the fluid density, p is pressure, x is distance along the pipe measured from the upstream, g is the acceleration of gravity, θ is the angle between the pipe and the x-axis, f is friction factor, D is the pipe diameter, and a is the wave speed. Equation (1) is the equation of motion and Equation (2) is the continuity equation for a compressible liquid in an elastic pipe.
These equations are complemented with a set of boundary conditions and require a numerical scheme to be solved. The method of characteristics (MOC) for the solution of the set of equations was utilized by the authors. Reservoirs, junctions, pumps, and valves are incorporated in the MOC as boundary conditions. Streeter & Wylie (1983) provide detailed derivation of the corresponding equations. The MOC is solved on a dense numerical grid in the spatial and temporal domains, defined by Δt and Δx, which are selected so that Δxa Δt. The optimization algorithm embeds a transient simulation model and requires one or several transient simulation runs each time the objective function is evaluated. Thus, the time step directly affects overall running times. At the same time, a grid step should be small enough to ensure numerical scheme stability. Balancing precision and computational time, the time step Δt = 2·10−2 s was chosen. It results in relative mistakes in transient simulation results of the order 10−5. Temporal discretization analysis of this transient simulation model for similar networks can be found in Skulovich et al. (2015). In this study, we consider control valves that serve to regulate the flow in the system by adjusting the valve opened area. Each type of valve, e.g., butterfly, gate, and needle valve, has a characteristic closing pattern, which controls the way the flow changes with respect to the change in the valve opened area. The change in the flow can be modeled using valve discharge coefficient Cv (Streeter & Wylie 1983): 
formula
3
where Qi,j is the flow at node i at time j, Cv is the discharge coefficient, Aj is the valve open area at time j, and Hi,j is the head at node i at time j. The value of the discharge coefficient Cv is determined for steady state and assumed constant (Streeter & Wylie 1983), while the changes are expressed by change in Aj.

Equation (3) is introduced as the controlled boundary condition at the location of the control valve and affects system response depending on the change in its cross-area. Using Equation (3), we assume that change in the flow is known for a given valve opening and we look for the optimal way to change the valve opening parameter, open area, with time. In the same manner, if another opening parameter defines change in the flow and can be controlled (e.g., stem travel distance), its change with time should be treated as a decision variable.

Examples for the change in the flow as a function of valve open area for different discharge coefficients Cv are shown in Figure 2. For example, for a valve with Cv = 0.269 an opening of only 20% of its area will result in 60% of its maximum capacity flow.
Figure 2

Valve characteristic curves.

Figure 2

Valve characteristic curves.

Optimization problem

To ensure minimal transient effect on the system, we examine how different valve closure strategies influence the resulting pressures, in other words, we are looking for the optimal way to change valve opening area with time to control the pressures in the system.

Given network parameters K (e.g., layout, pipes’ length, diameter, friction coefficient, demands at nodes), a set of equations describing system response to transients induced by operation of valves TS, a characteristic curve describing the area-flow valve response, and a control variable c, the objective of the optimization problem is to find the optimal control parameters such that the maximum pressure head in the system, z, is minimized. The problem is formulated as follows: 
formula
4

Valve characteristic curves

Three types of curves are examined as potential valve closure curves: quadratic, power, and piece-wise linear.

Quadratic formulation

Percent of valve opened area is defined as a second-order polynomial of the closure time: 
formula
5
where Aj is the valve open area at time tj, A0 is the initial valve open area, tcl is the total valve closure time and c is the curve coefficient.
With an assumption that at the beginning (t=tj=0 = 0) the valve is fully open (Aj=0 = A0) and at the end of closure time the valve is fully closed (Atcl = 0), the closure curve can be defined as a parabola with one unknown coefficient c. To ensure that open area fraction stays within (0,1) limits, coefficient c should be strictly limited by c ∈ [–1, 1]. Figure 3(a) shows the valve open area Aj as a function of control parameter for different values of c. Skulovich et al. (2014) used quadratic closure curve formulation with QN and GA optimization methods to minimize the transient resulted from valve closure in simple water distribution networks. The present work is a further development of the mentioned study, broadened with wider methodology and application examples.
Figure 3

Percent of valve opened area as a function of time for different values of curve coefficient c for (a) quadratic function, (b) power function, (c) piece-wise linear function.

Figure 3

Percent of valve opened area as a function of time for different values of curve coefficient c for (a) quadratic function, (b) power function, (c) piece-wise linear function.

Power formulation

The percent of valve opened area is defined as a power function: 
formula
6
A similar formulation can be found in Azoury et al. (1986). Example curves are given in Figure 3(b). This formulation can be viewed as the expansion of the first formulation with wider curve range available. To ensure that open area fraction stays within the required limits, coefficient c should be strictly positive: c > 0.

Piece-wise linear formulation

The percent of valve opened area is defined as a linear function for five regions with step 0.2 in closure time. Initially, this defines ten coefficients, but taking border conditions and curve continuity into account, only four decision variables, denoted as c1, c2, c3, and c4, are required.

If the five linear equations are denoted as: 
formula
7
then , , , and all other coefficients are derived from them.

Example curves are given in Figure 3(c). To ensure that the function remains descending, all the coefficients c1… c4 should be non-positive, for the function to remain within the defined rough limits, c1 >5, c2… c4 > 2.

Optimization methods

Two optimization approaches were considered in this study: genetic and QN algorithms, as described below.

GA

GA is a commonly used evolutionary algorithm for solving nonlinear nonsmooth optimization problems (Man et al. 1996). The algorithm mimics the process of natural evolution. A set of candidate solutions forms the population for the optimization problem. Each candidate solution in the population has a set of properties that can be mutated and altered. Starting from a certain initial population, the solutions iteratively change towards optimum solution improvement. This approach does not require certain restrictive conditions on the objective and the constraint functions (e.g., continuity and differentiability) which are not observed for transient states.

To solve the depicted valve closure problem, GA is applied directly using ga function in Matlab (MATLAB 2014) with decision variable c corresponding to the closing curve (one or four variables, depending on the chosen curve).

QN algorithm

QN methods are iterative and generally involve computations of objective function z and its gradient ∇z at each iteration. The basic form of the QN method is the following:

  1. Set number of iterations k=0. Make an initial guess ck. Initialize Hk = I, the n by n identity matrix, n is the number of decision variables.

    The matrix Hk is an approximation of the inverse of the Hessian at the current iterate. By approximation of the inverse of the Hessian we do not need to solve the linear system to find search direction, but only multiply this approximated matrix by the gradient (Bryan 2011).

  2. Compute objective function z(ck), the maximum pressure head in the system for given initial guess, and its numerical gradient ∇z(ck) and take the search direction hk as hk = −Hkz(ck). Numerical gradient is calculated using transient simulation to compute the response z for ck+ɛ and ckɛ. Then, numerical derivative (Davis & Polonsky 1972) 
    formula
    8
    where ci is the element of the vector ck, i∈ 1. n, n= 1 or n= 4 ɛ is taken equal to 10−4.
  3. Update ck = ck−1+hk, objective function, numerical gradient, and matrix Hk. In this study, the matrix update at each iteration from Hk−1 to Hk is the DFP (Davidson, Fletcher, and Powell) update. It takes Hk = Hk−1 + Uk, where 
    formula
    9
    where . This updates globally and superlinearly convergent for linear problems and locally and superlinearly convergent for nonlinear problems (Dennis & Moré 1977).
  4. Check convergence. The algorithm stops if the value of the objective function at the current step is different from the value of the objective function at the previous step for no more than ɛ = 10−4. If not, compute new search direction and go to step 3.

The QN method is the method for unconstrained optimization; therefore, decision variable bounds were not applied in this approach.

EXAMPLE APPLICATIONS

Three case studies were explored in this work through base runs and sensitivity analyses, as described below.

Case study 1

This example represents a simple WDS, comprising a reservoir with a constant head, one pipe, and a valve at the downstream end (Figure 4(a)). All elevations are zero. Head at the reservoir is 30 m, length of the pipe is 640 m, diameter is 900 mm, Darcy–Weisbach coefficient is 0.02, and wave speed is 1,000 m/s. Steady state velocity is 1.3 m/s.
Figure 4

Network layout: (a) case study 1, (b) case study 2.

Figure 4

Network layout: (a) case study 1, (b) case study 2.

Three valve closure curves and two optimization algorithms were applied to find the minimum of the maximum pressure head resulting from a 5-second valve closure at node 2.

Table 1 summarizes results for all six possible closure curve–optimization method configurations and provides the best maximum pressure head and average maximum pressure head with standard deviation from ten runs of each algorithm and each closing curve. For the QN method, initial guess was chosen randomly from the decision variables domain. The values of closure curves’ coefficients imply almost linear change in the valve opened area with time (Figure 5(a)). Quadratic and power curves are very close. Piece-wise linear curve suggests slightly faster closure at the beginning (c1 = −1.49 < −1) and slower closure at the end (c2… c4 >− 1).
Table 1

Case study 1: optimization results

Closure curve A=f(t) Number and domain of decision variables Optimal value of decision variables GA best Hmax (m) GA average Hmax ± stand. dev. (m) QN best Hmax(m) QN average Hmax± st. dev. (m) 
Quadratic 1, c ∈ (–1, 1) c= 0.127 53.05 53.05 ± 0 53.05 53.05 ± 0 
Power 1, c > 0 c= 1.041 53.71 53.71 ± 0 53.71 53.71 ± 0 
Piece-wise linear 4, c1 ∈ (–5, 0), c2–c4 ∈ (–2, 0) −1.49 −0.87 −0.89 −0.89 50.05 50.09 ± 0.04 50.05 87.21 ± 50.75 
Closure curve A=f(t) Number and domain of decision variables Optimal value of decision variables GA best Hmax (m) GA average Hmax ± stand. dev. (m) QN best Hmax(m) QN average Hmax± st. dev. (m) 
Quadratic 1, c ∈ (–1, 1) c= 0.127 53.05 53.05 ± 0 53.05 53.05 ± 0 
Power 1, c > 0 c= 1.041 53.71 53.71 ± 0 53.71 53.71 ± 0 
Piece-wise linear 4, c1 ∈ (–5, 0), c2–c4 ∈ (–2, 0) −1.49 −0.87 −0.89 −0.89 50.05 50.09 ± 0.04 50.05 87.21 ± 50.75 
Figure 5

GA and QN optimal valve closure curves: (a) case study 1, (b) case study 3.

Figure 5

GA and QN optimal valve closure curves: (a) case study 1, (b) case study 3.

Figure 6 shows 40-second transient simulation results for the pressure head at the valve. For example, a quadratic closure profile with coefficient c = −1 (high amplitude line) subjects the system to the repeated pressure rise to ∼90 m (three times higher than the steady state pressure of ∼30 m), along with severe pressure drop below cavitation head (lower than −10 m). The optimal solution obtained with GA for piece-wise linear closure curve results in slightly lower maximum pressure head than the maximum pressure head obtained with GA for quadratic and power closure curves. It can be noticed that the maximum pressure head occurs during the first 5 seconds, i.e., before the valve is completely closed. After the complete closure, the optimal results for all three closure curves are very close.
Figure 6

Case study 1: transient simulation results for different closure curves.

Figure 6

Case study 1: transient simulation results for different closure curves.

Figure 7 shows the convergence of the QN method as a function of the initial guess. It can be observed from the results that:
  • for quadratic and power closure curves, the QN method is numerically stable and converges to the same value from any initial guess within a variable domain (Figure 7(a) and 7(b));

  • for piece-wise linear closure curve, the QN method shows high sensitivity to the initial guess (Figure 7(c)).

Figure 7

QN convergence as a function of initial guess, for closure curves in (a) quadratic, (b) power, (c) piece-wise linear formulations.

Figure 7

QN convergence as a function of initial guess, for closure curves in (a) quadratic, (b) power, (c) piece-wise linear formulations.

Numerical sensitivity of the QN method can be explained by the convergence to local minimum for nonlinear problems such as transient. In addition, objective function in this study is a function comprising numerical values of the maximum pressure head, obtained from transient simulation. Its calculation makes sense only in a certain domain of the decision variables. During algorithm propagation, it is possible that current step size is larger than the decision variable domain, and further algorithm propagation is disturbed. QN method instability can be reduced by careful choice of the initial value for the approximation of the inverse Hessian matrix. Instead of the identity matrix, we suggest using H0 = 0.1·I, as any positive definite matrix can be used for the initial approximation (Bryan 2011). This change significantly widens the intervals of the initial guess that leads to the algorithm's convergence. For example, for case study 1, it ensures convergence for any initial guess from −3 to −0.4, whereas without the change this domain did not lead to convergence (Figure 7(c)).

An additional check was conducted to investigate whether a smaller number of parameters of the piece-wise closure curve leads to the same optimal result and better QN method convergence. One- and two-parameter formulations do not have enough flexibility to achieve the optimal maximum pressure head. Three-parameter formulation allows the optimal solution to be obtained, but QN algorithm convergence does not improve. Thus, four-parameter piece-wise linear closure curve without changes was chosen to be used in further case studies.

Sensitivity analysis

In the previous section, we presented optimization results for valve closure strategy that incudes closure curve combined with valve curve. In every case closure curve was combined with one and the same valve curve, valid for the data of the given example (Cv = 0.056). To check whether closure curve choice and optimization method choice are dependent on the valve curve, we repeat the same study for different discharge coefficients (obtained by change in the pipe diameter, keeping the steady state flow and head constant). The results are presented in Table 2. For each combination of the parameters, calculations were repeated five times, and average and standard deviation for the maximum pressure head are given. Smaller pipe diameter (higher discharge coefficient) corresponds to more severe transient event and results in higher maximum pressure head. This makes a more noticeable difference between the optimization results obtained with three closure curves. Power closure curve results are up to 7.9% higher and piece-wise linear closure curve results are up to 50.6% lower than the maximum pressure head value obtained with quadratic closure curve. Thus, power closure curve formulation does not provide any advantage despite a much wider range of allowable counters in comparison with quadratic curve. For both one-parametric closure curves GA and QN methods yield the same maximum pressure head value with very low deviation. As before, QN method combined with piece-wise linear closure curve is very sensitive to the initial guess, which explains the large standard deviation.

Table 2

Case study 1 sensitivity analysis: optimization results

Discharge coefficient Cv Pipe diameter (mm) Closure curve GA best Hmax (m) GA average Hmax± stand. dev. (m) QN best Hmax (m) QN average Hmax± st. dev. (m) 
0.42 500 Quadratic 193.49 193.49 ± 0a 193.49 193.49 ± 0a 
Power 201.41 201.41 ± 0a 201.41 201.41 ± 0a 
Piece-wise linear 95.58 96.03 ± 0.47 95.55 276 ± 196 
0.15 600 Quadratic 97.58 97.58 ± 0a 97.58 97.58 ± 0a 
Power 105.32 105.32 ± 0a 105.32 105.32 ± 0a 
Piece-wise linear 74.26 74.89 ± 1.07 96.14 202 ± 152 
0.099 700 Quadratic 73.29 73.29 ± 0 73.29 73.29 ± 0 
Power 76.93 76.93 ± 0 76.93 76.93 ± 0 
Piece-wise linear 62.73 63.0 ± 0.43 62.94 292 ± 175 
0.045 1,000 Quadratic 47.90 47.90 ± 0 47.90 47.90 ± 0 
Power 48.09 48.09 ± 0 48.09 48.09 ± 0 
Piece-wise linear 46.33 46.49 ± 0.31 46.39 145 ± 122 
0.019 1,400 Quadratic 38.56 38.56 ± 0a 38.56 38.56 ± 0a 
Power 38.61 38.61 ± 0a 38.61 38.61 ± 0a 
Piece-wise linear 38.58 38.64 ± 0.12 38.44 46.88 ± 17.4 
Discharge coefficient Cv Pipe diameter (mm) Closure curve GA best Hmax (m) GA average Hmax± stand. dev. (m) QN best Hmax (m) QN average Hmax± st. dev. (m) 
0.42 500 Quadratic 193.49 193.49 ± 0a 193.49 193.49 ± 0a 
Power 201.41 201.41 ± 0a 201.41 201.41 ± 0a 
Piece-wise linear 95.58 96.03 ± 0.47 95.55 276 ± 196 
0.15 600 Quadratic 97.58 97.58 ± 0a 97.58 97.58 ± 0a 
Power 105.32 105.32 ± 0a 105.32 105.32 ± 0a 
Piece-wise linear 74.26 74.89 ± 1.07 96.14 202 ± 152 
0.099 700 Quadratic 73.29 73.29 ± 0 73.29 73.29 ± 0 
Power 76.93 76.93 ± 0 76.93 76.93 ± 0 
Piece-wise linear 62.73 63.0 ± 0.43 62.94 292 ± 175 
0.045 1,000 Quadratic 47.90 47.90 ± 0 47.90 47.90 ± 0 
Power 48.09 48.09 ± 0 48.09 48.09 ± 0 
Piece-wise linear 46.33 46.49 ± 0.31 46.39 145 ± 122 
0.019 1,400 Quadratic 38.56 38.56 ± 0a 38.56 38.56 ± 0a 
Power 38.61 38.61 ± 0a 38.61 38.61 ± 0a 
Piece-wise linear 38.58 38.64 ± 0.12 38.44 46.88 ± 17.4 

aStandard deviation is ≤10−4.

Algorithms’ performance is compared using number of function evaluations, i.e., number of times the optimization algorithm refers to transient simulation routine, as a reference parameter.

GA requires 2,600–3,600 function evaluations when combined with quadratic closure curve, 2,900–3,900 function evaluations when combined with power closure curve, and 4,000–20,000 function evaluations when combined with piece-wise linear closure curve.

QN method requires 250–750 function evaluations when combined with quadratic closure curve, 100–1,100 function evaluations when combined with power closure curve, and 1,400–11,900 function evaluations when combined with piece-wise linear closure curve.

Thus, the QN method bears computational advantage in every case, but suffers from numerical sensitivity when coupled with piece-wise linear closure curve formulation. At the same time, piece-wise linear closure curve allows reduction of maximum pressure head more effectively than any other formulation, especially for severe transient case. In this case, GA might avoid local minima that cause a convergence issue for the QN method and, consequently, GA outperforms it.

Case study 2

The second case study was adopted from Streeter & Wylie (1967). The network comprises nine pipes, five junctions, a reservoir, three closed loops, and a valve located at the downstream end of the system, as shown in Figure 4(b). Head at the reservoir is constant and equal to Hres = 170 m. System data were slightly changed compared to the initial example to create more challenging conditions with cavitation risk for valve closure event. The data are listed in Table 3. The transient event is created by closing the valve located at node 7 within 2 seconds.

Table 3

Case study 2 data

Pipe Length (m) Diameter (mm) Friction factor Flow (steady state) (L/s) Wave speed (m/s) 
610 915 0.03 849.5 1,000 
915 760 0.028 409.6 1,000 
610 610 0.024 439.9 1,200 
550 460 0.02 176.5 900 
460 460 0.02 233.1 1,200 
490 460 0.025 112.0 1,000 
670 760 0.04 504.3 1,000 
460 610 0.03 345.2 900 
610 915 0.024 849.5 1,000 
Pipe Length (m) Diameter (mm) Friction factor Flow (steady state) (L/s) Wave speed (m/s) 
610 915 0.03 849.5 1,000 
915 760 0.028 409.6 1,000 
610 610 0.024 439.9 1,200 
550 460 0.02 176.5 900 
460 460 0.02 233.1 1,200 
490 460 0.025 112.0 1,000 
670 760 0.04 504.3 1,000 
460 610 0.03 345.2 900 
610 915 0.024 849.5 1,000 

Based on the results obtained for case study 1, power closure curve formulation was discarded as the least effective and only quadratic and piece-wise linear closure curves were explored.

Figure 8 shows 40-second transient simulation results for the pressure head at the valve. Instantaneous closure (noisy line) subjects the system to a repeated pressure rise up to 391.7 m and pressure drop far below cavitation head (to −61.9 m). The GA optimal solutions obtained with piece-wise linear closure curve and quadratic closure curve are very close. They result in maximum pressure head equal to 301.6 m and 302.4 m and corresponding minimum pressure head equal to 42.0 m and 48.4 m for piece-wise linear and quadratic closure curves, respectively. Optimal curve parameters are c= 0.297 (quadratic) and c1–c4 = (−1.38, −0.83, −0.97, −0.99) (piece-wise linear).
Figure 8

Case study 2: transient simulation results for different closure curves.

Figure 8

Case study 2: transient simulation results for different closure curves.

The QN method coupled with quadratic curve converges to similar results as attained by using GA results and is stable for any initial guess. The QN method requires 160–800 transient simulation runs whereas GA requires 2,600.

The QN method coupled with piece-wise linear closure curve results in even lower maximum pressure head of 278.3 m and corresponding higher minimum pressure head of 61.2 m. This result is not stable with respect to initial guess, but it is for the value of the initial guess larger than −0.5 for all decision variables: average maximum pressure head in five runs is 278.36 m with standard deviation only 0.07 m. In addition, this result is obtained within 1,000–2,800 transient simulation runs, whereas GA propagation requires 4,600–16,000 transient simulation runs.

Sensitivity analysis

The sensitivity of the optimization methods was explored by modifying case study 2 with the following changes:
  • Case study 2.1. Constant demand of 100 L/s is added at node 3 (Figure 9(a)).

  • Case study 2.2. Reservoir with constant head of 167 m is added at node 5 with new connecting pipe 10; all pipe parameters repeat the parameters of pipe 9 (Figure 9(b)).

  • Case study 2.3. New junction 8 is added in the middle of pipe 9, now dividing it into pipe 9 (new) and pipe 11; all parameters are equivalent to initial pipe 9. New pipe 10 connects node 5 and node 8; all pipe parameters repeat the parameters of pipe 8 (Figure 9(c)).

  • Case study 2.4. Change in diameter of pipe 9. Case 2.4.1: diameter is decreased to 0.61 m; case 2.4.2: diameter is increased to 1.3 m (Figure 9(d)).

Figure 9

Case study 2: sensitivity analysis. Network layouts for cases (a) 2.1, (b) 2.2, (c) 2.3, (d) 2.4.

Figure 9

Case study 2: sensitivity analysis. Network layouts for cases (a) 2.1, (b) 2.2, (c) 2.3, (d) 2.4.

GA and QN optimization methods were used to find optimal quadratic and piece-wise linear closure curves. The obtained maximum and minimum pressure head during 20 second transient simulation for every optimal closure curve in comparison with instantaneous closure is given in Table 4. For each combination of parameters, the algorithm was executed five times for GA and five times from different initial guesses for QN. The obtained results are either stable with standard deviation lower than 10−4, or highly unstable (for piece-wise linear) with convergence to optimum only from narrow interval of the initial guess.

Table 4

Optimization results for cases 2–2.4.2

Case study Instantaneous closure
 
GA: optimal closure (quadratic)
 
QN: optimal closure (quadratic)
 
GA: optimal closure (piece-wise linear)
 
QN: optimal closure(piece-wise linear)
 
Hmax (m) Hmin(m) Hmax(m) Hmin (m) Hmax (m) Hmin(m) Hmax (m) Hmin (m) Hmax (m) Hmin (m) 
391.7 −61.9 302.4 48.4 302.4 48.4 301.6 42.0 278.4a 61.2 
2.1 390.9 −62.9 301.7 48.1 301.7 48.1 301.5 42.2 300.5b 41.8 
2.2 339.3 3.7 286.1 61.9 286.1 61.9 279.5 42.6 279.4b,c 40.3 
2.3 389.7 −29.9 283.9 57.9 283.9 57.9 284.0 57.7 283.4c 54.9 
2.4.1 571.9 −206.6 376.2 −51.6 376.2c −51.6 339.1 −12.1 337.9c −9.7 
2.4.2 340.1 −46.6 279.1 59.4 279.1 59.4 279.0 59.5 278.3b 61.2 
Case study Instantaneous closure
 
GA: optimal closure (quadratic)
 
QN: optimal closure (quadratic)
 
GA: optimal closure (piece-wise linear)
 
QN: optimal closure(piece-wise linear)
 
Hmax (m) Hmin(m) Hmax(m) Hmin (m) Hmax (m) Hmin(m) Hmax (m) Hmin (m) Hmax (m) Hmin (m) 
391.7 −61.9 302.4 48.4 302.4 48.4 301.6 42.0 278.4a 61.2 
2.1 390.9 −62.9 301.7 48.1 301.7 48.1 301.5 42.2 300.5b 41.8 
2.2 339.3 3.7 286.1 61.9 286.1 61.9 279.5 42.6 279.4b,c 40.3 
2.3 389.7 −29.9 283.9 57.9 283.9 57.9 284.0 57.7 283.4c 54.9 
2.4.1 571.9 −206.6 376.2 −51.6 376.2c −51.6 339.1 −12.1 337.9c −9.7 
2.4.2 340.1 −46.6 279.1 59.4 279.1 59.4 279.0 59.5 278.3b 61.2 

aThe QN algorithm converges in a wide interval of the initial guess, for example [−0.5, 1].

bThe QN algorithm converges only in a narrow interval of the initial guess, for example [−0.19, −0.17].

cAt a certain point the valve opens back slightly.

Case 2.1 is similar to the initial case study 2 both in instantaneous closure results and optimal results. This is to be expected since adding extra loading to the network does not influence transient behavior. Case 2.2 is the mildest transient event with lower impact on the system than the base case 2, because additional reservoir acts like a surge tank mitigating the transient. Case 2.3 is a moderate severity case. In both cases 2.2 and 2.3, new reservoir and pipe (pipe 10) provide additional means for wave energy to dissipate and, therefore in these cases, the maximum pressure head is lower and the minimum pressure head is slightly higher than in the initial case study 2. Significant diameter decrease (case 2.4.1) results in very severe transient event. Consequently, for this event, in particular, the optimization results for quadratic and piece-wise linear closure curves vary the most. Although optimal closure according to the quadratic curve significantly reduces the maximum pressure head (34% reduction in comparison with instantaneous closure) and significantly increases the minimum pressure head (75% increase), the minimum pressure head is still very low and far below cavitation head. At the same time, closure according to the optimal piece-wise linear closure curve provides significantly better protection and almost eliminates cavitation. Case 2.4.2 is a mild transient event and the difference between two closing curves is least noticeable.

On average, both quadratic and piece-wise linear closure curves obtained with GA reduce the maximum pressure head by 23.4 and 26.3%, respectively; in other words, the difference between the two closing curves is minor. Slight difference in the results can be noticed for case 2.2. Significant difference in the results from two types of curves is evident only for severe transient case 2.4.1.

GA optimal closure curves for the denoted cases are presented in Figure 10. Interestingly, for most cases, the tendencies are similar for the two types of curves. The most evident difference in the shape is for case 2.4.1 – the most severe transient event. Here, the flexible structure of the piece-wise linear formulation allows change in the open area in an acute way that, in the end, results in lower maximum pressure head and higher minimum pressure head.
Figure 10

GA optimal closure curves (a) quadratic and (b) piece-wise linear for case study 2.

Figure 10

GA optimal closure curves (a) quadratic and (b) piece-wise linear for case study 2.

Comparing GA and QN methods, the following observations can be made:

  1. For quadratic (one-parameter) closure curve, QN results coincide with GA results. For this closure curve the QN method is stable and converges for any initial guess, except for severe case 2.4.1. For case 2.4.1, only a narrow interval of the initial guess led to optimal results for the QN method.

  2. For piece-wise linear (four-parameter) closure curve, QN results are slightly better than GA results, but the QN method suffers from instability. For relatively mild cases (2, 2.2 and 2.4.2) QN algorithm converges if initial guess is larger or smaller than certain value within the parameters’ domain. For other cases, only certain values of the initial guess led to the optimal solution. For example, for case 2.1, the initial guess should be within the interval (−0.19, −0.17) for all the parameters.

  3. By definition, the QN method is the method for unconstrained optimization. Therefore, the values of the decision variables are not limited in this method. It led to interesting results in case study 2.2. The solution obtained with GA suggests a pause in closure (Figure 10(b), the lowest line), since the decision variables are defined as non-positive (to define valve closure as non-returnable). QN lacks that definition; therefore, the optimal curve obtained with this method includes partial reopening of the valve before complete closure.

  4. From a computational aspect, the previously found tendency holds: QN requires two to ten times less transient simulation runs than GA.

Case study 3

The network comprises 11 pipes arranged in two loops as shown in Figure 11, two reservoirs at nodes 1 and 10, a pump at node 2, and a valve at node 9. Reservoir head at node 1 is constant and equal to 30 m, reservoir head at node 10 is constant and equal to 59.3 m. Two constant demands of 56.6 L/s and 28.3 L/s are defined at nodes 5 and 6, respectively. All elevations are zero. Wave speed is 1,200 m/s for all pipes. Pump curve is defined by three nodes: for zero flow the change in head is 30 m, for 84.9 L/s flow the change in head is 24.4 m, and for 169.9 L/s flow the change in head is 15 m. Additional network data are given in Table 5.
Figure 11

Network layout: case study 3.

Figure 11

Network layout: case study 3.

Table 5

Case study 3 data

Pipe Length (m) Diameter (mm) Friction factor 
30 300 0.0379 
60 300 0.0379 
45 250 0.04 
45 250 0.04 
45 200 0.0437 
45 200 0.0437 
60 200 0.0437 
60 200 0.0437 
60 200 0.0437 
10 60 200 0.0437 
11 30 200 0.0437 
Pipe Length (m) Diameter (mm) Friction factor 
30 300 0.0379 
60 300 0.0379 
45 250 0.04 
45 250 0.04 
45 200 0.0437 
45 200 0.0437 
60 200 0.0437 
60 200 0.0437 
60 200 0.0437 
10 60 200 0.0437 
11 30 200 0.0437 

During the steady state, the maximum pressure head is equal to 62.1 m and the minimum pressure head is 39.6 m.

Figure 12(a) shows 30-second transient simulation results for the pressure head at the valve after 1-second valve closure. Inappropriate valve closure subjects the system to a pressure up to 132.9 m and down to −9.9 m, introducing cavitation risk.
Figure 12

Case study 3: transient simulation results for inappropriate and GA and QN optimal (a) valve closure and (b) pump shutdown.

Figure 12

Case study 3: transient simulation results for inappropriate and GA and QN optimal (a) valve closure and (b) pump shutdown.

GA optimization results in the maximum pressure head equal to 95.2 m and 76.3 m and corresponding minimum pressure head equal to 34.6 m and 37.4 m for quadratic and piece-wise linear closure curves, respectively. Optimal curve parameters are c= 0.94 (quadratic) and c1–c4 = (−3.89, −0.64, −0.07, −0.16) (piece-wise linear), and the curves are shown in Figure 5(b). QN results for quadratic curve repeat the GA results; the method is stable for initial guess lower than 0.1 and requires less transient simulation runs than GA. QN results for piece-wise linear closure curve are 77.0 m maximum pressure head and 35.4 m minimum pressure head; in other words, QN results for this case are no better than GA results. In addition, the method shows high instability and converges to the denoted result only from certain values of the initial guess.

Pump shutdown

In this case, we apply closure curves and optimization methodology to find an optimal way to shut down the pump at node 2 within 4 seconds so that it results in the minimal maximum pressure head. For a variable speed pump, its characteristic curve as a function of the speed can be formulated as (Streeter & Wylie 1983): 
formula
10
where N is the pump speed (rounds per minute) and C1C3 are the pump curve coefficients, H is the hydraulic grade line, and Q is the discharge.

Then we reformulate closure curves from previous case studies as change in speed curves, meaning quadratic and piece-wise curves now describe the change in pump speed N within pump stop time. Standard transient modeling methodology is applied to find all other parameters.

During controlled pump shutdown, a linear speed change is often assumed for the speed variation. Under this assumption, a 4-second pump shutdown subjects the system to a pressure head up to 102.6 m and pressure drop to −25.3 m, introducing cavitation. GA optimization results are:

  • for quadratic speed change: maximum pressure head 79.6 m, minimum pressure head 3.8 m, optimal curve coefficient c = 1;

  • for piece-wise linear speed change: maximum pressure head 69.3 m, minimum pressure head 24.3 m, optimal curve coefficients c1 = −4.97, c2c4 = 0.

QN results are slightly better:

  • for quadratic change in the speed: 74.1 m and 13.6 m minimum and maximum pressure head, respectively, optimal curve coefficient c = 2.02;

  • for piece-wise linear change in the speed: 68.5 m and 24.1 m minimum and maximum pressure head, respectively, optimal curve coefficients c1 = −6.42, c2 = 3.53, c3 = − 2.96; c4 = 1.12.

Transient simulation results for this case and QN optimal results are presented in Figure 12(b). In this case, the QN method is stable with respect to initial guess. However, since the method is unconstrained, the resulting coefficients for speed change curves imply secondary speeding up of the pump after the shutdown process is started before the complete pump stop (Figure 12(b), oscillations within first 4 seconds).

The obtained results imply that appropriate change in speed allows avoiding cavitation and negative head in the system and reduces the maximum pressure head significantly. Hence, without change in the pump shutdown time, negative effects of the transient can be eliminated by choosing the appropriate way to change pump speed only.

CONCLUSIONS

Rapid changes in operation of valves and pumps can result in significant pressure head oscillations above the maximum design head and below cavitation head. When increasing operational time of control devices is unavailable, the negative effects can be mitigated if proper closure strategy is applied. Three types of valve closure curves were explored: power, quadratic, and piece-wise linear function. Among proposed closure curves, piece-wise linear formulation is preferable as the one resulting in the minimal maximum pressure head, especially for severe transient events. Quadratic formulation is the more favored choice than power formulation; it results in sufficiently good solution and it is computationally advantageous in comparison with piece-wise linear formulation.

QN optimization method requires less transient simulation runs and for some cases results in better solutions than GA. The QN method is sensitive to the selection of initial guess, particularly for severe transient cases and for multi-parameter closure curves. The sensitivity can be partially reduced by careful selection of the initial approximation of the inverse of the Hessian matrix.

The obtained values of the curve coefficients vary over a wide range within the parameter's domain (for example, from 0.16 to 0.94 for quadratic closure curve). Therefore, the use of the proper optimization method is more effective than general rules of thumb. At the same time, all the curves imply faster closure at the beginning and slower closure at the end, in compliance with known closure pattern (Wood et al. 2005).

The proposed framework was shown to be applicable for pump shutdown scheduling.

The results are demonstrated for networks with a relatively small number of nodes. Since optimization methods rely on simulation of the transient conditions, which is known to be computationally extensive, scaling up to a large system is a challenging task. Thus, further work should be aimed at methods that require less function evaluations.

Future work will focus on: (1) improving the QN algorithm to reduce its instability by investigating sensitivity to the initial guess, initial Hessian matrix choice, step size choice; and (2) considering the use of the optimization framework for systems with multiple valves and other operational scheduling (such as pump start-up, fire hydrant flushing) in WDSs subjected to transient.

ACKNOWLEDGEMENTS

This study was supported by the Technion Grand Water Research Institute and by the joint Israeli Office of the Chief Scientist (OCS) Ministry of Science, Technology and Space (MOST), and by the German Federal Ministry of Education and Research (BMBF), under project no. 02WA1298.

REFERENCES

REFERENCES
Azoury
P. H.
Baasiri
M.
Najm
H.
1986
Effect of valve-closure schedule on water hammer
.
J. Hydraul. Eng.
112
(
10
),
890
903
.
Bazargan-Lari
M. R.
Kerachian
R.
Afshar
H.
Bashi-Azghadi
S. N.
2013
Developing an optimal valve closing rule curve for real-time pressure control in pipes
.
J. Mech. Sci. Technol.
27
(
1
),
215
225
.
Bergant
A.
Simpson
A. R.
Tijsseling
A. S.
2006
Water hammer with column separation: a historical review
.
J. Fluid. Struct.
22
(
2
),
135
171
.
Boulos
P. F.
Karney
B. W.
Wood
D. J.
Lingireddy
S.
2005
Hydraulic transient guidelines for protecting water distribution systems
.
J. AWWA
97
(
5
),
111
124
.
Bryan
K.
2011
.
Chaudhry
M.
1979
Applied Hydraulic Transients.
Van Nostrand Reinhold
,
New York, NY
,
USA
.
Davis
P. J.
Polonsky
I.
1972
Numerical Interpolation, Differentiation, and Integration. Handbook of Mathematical Functions
.
Dover Publications, Inc.
,
New York, NY
,
USA
.
Dennis
J. E.
Jr
Moré
J. J.
1977
Quasi-Newton methods, motivation and theory
.
SIAM Rev.
19
(
1
),
46
89
.
Goldberg
D. E.
Karr
C. L.
1987
Quick stroking: design of time-optimal valve motion
.
J. Hydraul. Eng.
113
(
6
),
780
795
.
Larock
B. E.
Jeppson
R. W.
Watters
G. Z.
2000
Hydraulics of Pipeline Systems.
CRC Press
,
Boca Raton, FL
,
USA
.
Man
K. F.
Tang
K. S.
Kwong
S.
1996
Genetic algorithms: concepts and applications
.
IEEE Trans. Indust. Electr.
43
(
5
),
519
534
.
MATLAB
2014
MATLAB and Optimization Toolbox, Release R2014a
.
The MathWorks, Inc.
,
Natick, MA
,
USA
.
Merritt
H. E.
1967
Hydraulic Control Systems
.
John Wiley & Sons
,
New York, NY
,
USA
.
Skulovich
O.
Perelman
L.
Ostfeld
A.
2014
Modeling and optimizing hydraulic transients in water distribution systems
.
Procedia Eng.
70
,
1558
1565
.
Skulovich
O.
Bent
R.
Judi
D.
Perelman Sela
L.
Ostfeld
A.
2015
Piece-wise mixed integer programming for optimal sizing of surge control devices in water distribution systems
.
Water Resour. Res.
51
(
6
),
4391
4408
.
Streeter
V. L.
1966
Valve stroking for complex piping systems
.
J. Hydraul. Div.
93
(
3
),
81
98
.
Streeter
V. L.
Wylie
E. B.
1967
Hydraulic Transients.
McGraw-Hill
,
New York
,
USA
.
Streeter
V. L.
Wylie
E. B.
1983
Fluid Mechanics: First SI Metric Edition.
McGraw Hill
,
New York
,
USA
.
Wood
D. J.
Boulos
P. F.
Lingireddy
S.
2005
Pressure Wave Analysis of Transient Flow in Pipe Distribution Systems
,
1st edn
.
MWH Soft
,
Pasadena, CA
,
USA
.