Optimal closure of system actuators for transient control: an analytical approach

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.

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.

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.

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

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₁, c₂, c₃, and c₄, are required.

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

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

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 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 c_k. Initialize H_k = I, the n by n identity matrix, n is the number of decision variables.

   The matrix H_k 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(c_k), the maximum pressure head in the system for given initial guess, and its numerical gradient ∇z(c_k) and take the search direction h_k as h_k = −H_k∇z(c_k). Numerical gradient is calculated using transient simulation to compute the response z for c_k+ɛ and c_k − ɛ. Then, numerical derivative (Davis & Polonsky 1972)

   

   8

   where c_i is the element of the vector c_k, i∈ 1. n, n= 1 or n= 4 ɛ is taken equal to 10⁻⁴.
3. Update c_k = c_k−1+h_k, objective function, numerical gradient, and matrix H_k. In this study, the matrix update at each iteration from H_k−1 to H_k is the DFP (Davidson, Fletcher, and Powell) update. It takes H_k = H_k−1 + U_k, where

   

   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⁻⁴. 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.

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

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.

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

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.

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.

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.

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.

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.

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⁻⁴, or highly unstable (for piece-wise linear) with convergence to optimum only from narrow interval of the initial guess.

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:

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.

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₁–c₄ = (−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.

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₁ = −4.97, c₂… c₄ = 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₁ = −6.42, c₂ = 3.53, c₃ = − 2.96; c₄ = 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.

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.

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.