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

### Transient modeling

*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: 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.

*t*and Δ

*x*, which are selected so that Δ

*x*=±

*a*Δ

*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

*C*(Streeter & Wylie 1983): where

_{v}*Q*is the flow at node

_{i,j}*i*at time

*j*,

*C*is the discharge coefficient,

_{v}*A*is the valve open area at time

_{j}*j*, and

*H*is the head at node

_{i,j}*i*at time

*j*. The value of the discharge coefficient

*C*is determined for steady state and assumed constant (Streeter & Wylie 1983), while the changes are expressed by change in

_{v}*A*.

_{j}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.

*C*are shown in Figure 2. For example, for a valve with

_{v}*C*= 0.269 an opening of only 20% of its area will result in 60% of its maximum capacity flow.

_{v}### 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.

*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:

### Valve characteristic curves

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

#### Quadratic formulation

*A*is the valve open area at time

_{j}*t*,

_{j}*A*

_{0}is the initial valve open area,

*t*is the total valve closure time and

_{cl}*c*is the curve coefficient.

*t*

*=*

*t*

_{j=0}= 0) the valve is fully open (

*A*

_{j}_{=0}=

*A*

_{0}) and at the end of closure time the valve is fully closed (

*A*

_{tcl}= 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

*A*as a function of control parameter for different values of

_{j}*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.

#### Power formulation

*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 *c*_{1}, *c*_{2}, *c*_{3}, and *c*_{4}, are required.

Example curves are given in Figure 3(c). To ensure that the function remains descending, all the coefficients *c*_{1}*… c*_{4} should be non-positive, for the function to remain within the defined rough limits, *c*_{1} >*–*5, *c*_{2}*… c*_{4} > *–*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:

*Set number of iterations k**=**0.*Make an initial guess*c*. Initialize H_{k}= I, the_{k}*n*by*n*identity matrix,*n*is the number of decision variables.The matrix H

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)._{k}*Compute objective function z(c*, the maximum pressure head in the system for given initial guess, and its numerical gradient ∇_{k})*z*(*c*) and take the search direction_{k}*h*as_{k}*h*= −H_{k}∇_{k}*z*(*c*). Numerical gradient is calculated using transient simulation to compute the response_{k}*z*for*c*_{k}*+**ɛ*and*c*−_{k}*ɛ*. Then, numerical derivative (Davis & Polonsky 1972) where*c*is the element of the vector_{i}*c*,_{k}*i*∈ 1.*n*,*n**=*1 or*n**=*4*ɛ*is taken equal to 10^{−4}.*Update c*_{k}= c_{k−1}*+**h*, numerical gradient, and matrix H_{k}, objective functionIn this study, the matrix update at each iteration from H_{k}._{k}_{−1}to His the DFP (Davidson, Fletcher, and Powell) update. It takes H_{k}= H_{k}_{k}_{−1}+*U*, where where . This updates globally and superlinearly convergent for linear problems and locally and superlinearly convergent for nonlinear problems (Dennis & Moré 1977)._{k}*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

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.

Closure curve A=f(t) | Number and domain of decision variables | Optimal value of decision variables | GA best H (m) _{max} | GA average H ± stand. dev. (m) _{max} | QN best H(m) _{max} | QN average H± st. dev. (m) _{max} |
---|---|---|---|---|---|---|

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, c_{1} ∈ (–5, 0), c_{2}–c_{4} ∈ (–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 H (m) _{max} | GA average H ± stand. dev. (m) _{max} | QN best H(m) _{max} | QN average H± st. dev. (m) _{max} |
---|---|---|---|---|---|---|

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, c_{1} ∈ (–5, 0), c_{2}–c_{4} ∈ (–2, 0) | −1.49 −0.87 −0.89 −0.89 | 50.05 | 50.09 ± 0.04 | 50.05 | 87.21 ± 50.75 |

*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.

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)).

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 H_{0} = 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 (*C _{v}* = 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.

Discharge coefficient Cv | Pipe diameter (mm) | Closure curve | GA best H (m) _{max} | GA average H± stand. dev. (m) _{max} | QN best H (m) _{max} | QN average H± st. dev. (m) _{max} |
---|---|---|---|---|---|---|

0.42 | 500 | Quadratic | 193.49 | 193.49 ± 0^{a} | 193.49 | 193.49 ± 0^{a} |

Power | 201.41 | 201.41 ± 0^{a} | 201.41 | 201.41 ± 0^{a} | ||

Piece-wise linear | 95.58 | 96.03 ± 0.47 | 95.55 | 276 ± 196 | ||

0.15 | 600 | Quadratic | 97.58 | 97.58 ± 0^{a} | 97.58 | 97.58 ± 0^{a} |

Power | 105.32 | 105.32 ± 0^{a} | 105.32 | 105.32 ± 0^{a} | ||

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 ± 0^{a} | 38.56 | 38.56 ± 0^{a} |

Power | 38.61 | 38.61 ± 0^{a} | 38.61 | 38.61 ± 0^{a} | ||

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 H (m) _{max} | GA average H± stand. dev. (m) _{max} | QN best H (m) _{max} | QN average H± st. dev. (m) _{max} |
---|---|---|---|---|---|---|

0.42 | 500 | Quadratic | 193.49 | 193.49 ± 0^{a} | 193.49 | 193.49 ± 0^{a} |

Power | 201.41 | 201.41 ± 0^{a} | 201.41 | 201.41 ± 0^{a} | ||

Piece-wise linear | 95.58 | 96.03 ± 0.47 | 95.55 | 276 ± 196 | ||

0.15 | 600 | Quadratic | 97.58 | 97.58 ± 0^{a} | 97.58 | 97.58 ± 0^{a} |

Power | 105.32 | 105.32 ± 0^{a} | 105.32 | 105.32 ± 0^{a} | ||

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 ± 0^{a} | 38.56 | 38.56 ± 0^{a} |

Power | 38.61 | 38.61 ± 0^{a} | 38.61 | 38.61 ± 0^{a} | ||

Piece-wise linear | 38.58 | 38.64 ± 0.12 | 38.44 | 46.88 ± 17.4 |

^{a}Standard 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 *H _{res}* = 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.

Pipe | Length (m) | Diameter (mm) | Friction factor | Flow (steady state) (L/s) | Wave speed (m/s) |
---|---|---|---|---|---|

1 | 610 | 915 | 0.03 | 849.5 | 1,000 |

2 | 915 | 760 | 0.028 | 409.6 | 1,000 |

3 | 610 | 610 | 0.024 | 439.9 | 1,200 |

4 | 550 | 460 | 0.02 | 176.5 | 900 |

5 | 460 | 460 | 0.02 | 233.1 | 1,200 |

6 | 490 | 460 | 0.025 | 112.0 | 1,000 |

7 | 670 | 760 | 0.04 | 504.3 | 1,000 |

8 | 460 | 610 | 0.03 | 345.2 | 900 |

9 | 610 | 915 | 0.024 | 849.5 | 1,000 |

Pipe | Length (m) | Diameter (mm) | Friction factor | Flow (steady state) (L/s) | Wave speed (m/s) |
---|---|---|---|---|---|

1 | 610 | 915 | 0.03 | 849.5 | 1,000 |

2 | 915 | 760 | 0.028 | 409.6 | 1,000 |

3 | 610 | 610 | 0.024 | 439.9 | 1,200 |

4 | 550 | 460 | 0.02 | 176.5 | 900 |

5 | 460 | 460 | 0.02 | 233.1 | 1,200 |

6 | 490 | 460 | 0.025 | 112.0 | 1,000 |

7 | 670 | 760 | 0.04 | 504.3 | 1,000 |

8 | 460 | 610 | 0.03 | 345.2 | 900 |

9 | 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.

*c*

*=*0.297 (quadratic) and

*c*

_{1}

*–c*

_{4}= (−1.38, −0.83, −0.97, −0.99) (piece-wise linear).

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

*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)).

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.

Case study | Instantaneous closure | GA: optimal closure (quadratic) | QN: optimal closure (quadratic) | GA: optimal closure (piece-wise linear) | QN: optimal closure(piece-wise linear) | |||||
---|---|---|---|---|---|---|---|---|---|---|

H (m) _{max} | H(m) _{min} | H(m) _{max} | H (m) _{min} | H (m) _{max} | H(m) _{min} | H (m) _{max} | H (m) _{min} | H (m) _{max} | H (m) _{min} | |

2 | 391.7 | −61.9 | 302.4 | 48.4 | 302.4 | 48.4 | 301.6 | 42.0 | 278.4^{a} | 61.2 |

2.1 | 390.9 | −62.9 | 301.7 | 48.1 | 301.7 | 48.1 | 301.5 | 42.2 | 300.5^{b} | 41.8 |

2.2 | 339.3 | 3.7 | 286.1 | 61.9 | 286.1 | 61.9 | 279.5 | 42.6 | 279.4^{b,c} | 40.3 |

2.3 | 389.7 | −29.9 | 283.9 | 57.9 | 283.9 | 57.9 | 284.0 | 57.7 | 283.4^{c} | 54.9 |

2.4.1 | 571.9 | −206.6 | 376.2 | −51.6 | 376.2^{c} | −51.6 | 339.1 | −12.1 | 337.9^{c} | −9.7 |

2.4.2 | 340.1 | −46.6 | 279.1 | 59.4 | 279.1 | 59.4 | 279.0 | 59.5 | 278.3^{b} | 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) | |||||
---|---|---|---|---|---|---|---|---|---|---|

H (m) _{max} | H(m) _{min} | H(m) _{max} | H (m) _{min} | H (m) _{max} | H(m) _{min} | H (m) _{max} | H (m) _{min} | H (m) _{max} | H (m) _{min} | |

2 | 391.7 | −61.9 | 302.4 | 48.4 | 302.4 | 48.4 | 301.6 | 42.0 | 278.4^{a} | 61.2 |

2.1 | 390.9 | −62.9 | 301.7 | 48.1 | 301.7 | 48.1 | 301.5 | 42.2 | 300.5^{b} | 41.8 |

2.2 | 339.3 | 3.7 | 286.1 | 61.9 | 286.1 | 61.9 | 279.5 | 42.6 | 279.4^{b,c} | 40.3 |

2.3 | 389.7 | −29.9 | 283.9 | 57.9 | 283.9 | 57.9 | 284.0 | 57.7 | 283.4^{c} | 54.9 |

2.4.1 | 571.9 | −206.6 | 376.2 | −51.6 | 376.2^{c} | −51.6 | 339.1 | −12.1 | 337.9^{c} | −9.7 |

2.4.2 | 340.1 | −46.6 | 279.1 | 59.4 | 279.1 | 59.4 | 279.0 | 59.5 | 278.3^{b} | 61.2 |

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

^{b}The QN algorithm converges only in a narrow interval of the initial guess, for example [−0.19, −0.17].

^{c}At 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.

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

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.

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.

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.

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

### Case study 3

Pipe | Length (m) | Diameter (mm) | Friction factor |
---|---|---|---|

1 | 30 | 300 | 0.0379 |

2 | 60 | 300 | 0.0379 |

3 | 45 | 250 | 0.04 |

4 | 45 | 250 | 0.04 |

5 | 45 | 200 | 0.0437 |

6 | 45 | 200 | 0.0437 |

7 | 60 | 200 | 0.0437 |

8 | 60 | 200 | 0.0437 |

9 | 60 | 200 | 0.0437 |

10 | 60 | 200 | 0.0437 |

11 | 30 | 200 | 0.0437 |

Pipe | Length (m) | Diameter (mm) | Friction factor |
---|---|---|---|

1 | 30 | 300 | 0.0379 |

2 | 60 | 300 | 0.0379 |

3 | 45 | 250 | 0.04 |

4 | 45 | 250 | 0.04 |

5 | 45 | 200 | 0.0437 |

6 | 45 | 200 | 0.0437 |

7 | 60 | 200 | 0.0437 |

8 | 60 | 200 | 0.0437 |

9 | 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.

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 *c*_{1}*–c*_{4} = (−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

*N*is the pump speed (rounds per minute) and

*C*

_{1}–

*C*

_{3}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

*c*_{1}= −4.97,*c*_{2}…*c*_{4}= 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

*c*_{1}= −6.42,*c*_{2}= 3.53,*c*_{3}= − 2.96;*c*_{4}= 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.