## Abstract

Transients are used as a diagnostic tool for pressurized pipes due to their capabilities to acquire and transmit information about system status. To interpret such information in complex systems, models obtained by the integration of the governing equations in the frequency domain are used, since they are computationally efficient and reliable. These models do not require a regular spatial grid for integration and introduce the pipe lengths as parameters. In this paper a calibration procedure based on this particular feature is introduced, to determine a basic item of information that can be lost about the water distribution system, that is the length of the pipes. The network admittance matrix method is implemented and its numerical efficiency allows the investigation of the calibration optimization function on a regular grid in the parameter space. The calibration is tested using in series a genetic and the Nelder–Mead algorithms, considering both elastic and viscoelastic material pipes. The results of the presented numerical investigation allow some insights into the existence of the solution and into the shape of the optimization function.

## INTRODUCTION

Transients are a useful tool for the diagnosis of pressurized pipe systems. Starting from the pioneering works of the 1990s (Jonsson & Larson 1992; Liggett & Chen 1994; Liou & Tian 1995; Brunone 1999), different transient-test-based techniques have been developed and used so far. The proposed categorizations of these techniques consider the domain used for the integration of the transient governing equations (i.e. space and time or frequency domain), the inverse problem position (i.e. only the first wave arrival time or a general inverse transient analysis), and the detected anomalies (e.g. leaks and blockages).

Since techniques following different approaches are still used and proposed, it is evident that no single approach prevails over the others for all the possible applications.

In recent decades many authors have used models derived by the integration of the governing equations in the frequency domain, or frequency domain models (e.g. Ferrante & Brunone 2003; Lee *et al.* 2007; Sattar *et al.* 2008; Gong *et al.* 2014; Zecchin *et al.* 2014; Kim 2015; Ferrante *et al.* 2016; Duan 2017). These models have the main drawback of requiring a linearization of the governing equations. If the linearization error is negligible with respect to the other introduced approximations (Lee & Vítkovský 2010; Capponi *et al.* 2017b), they can be used with a high computational efficiency to simulate transients in pressurized pipe systems without the need of a fixed space grid. This advantage is fundamental for inverse analysis in complex systems, where computational burden and grid flexibility play a crucial role. Furthermore, with respect to classical time domain models, such as those based on the method of characteristics, the frequency domain integration introduces in the equations the length of the pipes as a parameter (Zecchin *et al.* 2013). For this reason, using the frequency domain models, lengths can be directly calibrated by means of an inverse analysis.

In the following, based on this feature of frequency domain integration, a pipe length determination problem in a network is considered. In practical applications it corresponds to the problem of buried and forgotten pipes, when the geometry of the system is not completely known because it has been forgotten, or when illegal connections of unknown length have been added.

With this aim, a branched pipe system is considered and a frequency domain model based on the network admittance matrix method (NAMM) is used. Such a system is relatively simple, but the introduction of a junction (Wood & Chao 1971) allows a step to be taken towards more complex systems and to test diagnosis methods on branched systems, which have not been widely analyzed in the literature yet. There are, in fact, few contributions in the field of transients in Y-pipe systems (Brunone *et al.* 2008; Ferrante *et al.* 2009; Meniconi *et al.* 2011; Evangelista *et al.* 2015; Capponi *et al.* 2017c).

The reduced computational burden of NAMM allows the evaluation on a regular grid of the optimization function used in the calibration procedure. As a further test, once the optimization function is known, a calibration procedure is applied combining in series a genetic and a non-linear optimization algorithm (NOA). Failures and successes of such calibration can be related to the main features of the known shape of the optimization function. Two numerical case studies, comprising elastic and viscoelastic pipe material systems, are considered.

The main contributions of this paper are the straightforward introduction of the lengths of the pipe as unknowns in a calibration procedure based on transient test data and the investigation of the optimization function by direct scrutiny, for a branched system. Moreover some interesting remarks are given about the influence of pipe material rheology on calibration performance and results.

## THE NAMM

*et al.*2009, 2010). It makes use of graph theoretic concepts to organize the 1-D water hammer governing equations, perturbed, linearized and integrated in the frequency domain, in an admittance matrix form: In Equation (1) the vectors

**and**

*Ψ***contain the Fourier transform of the variations of the pressure head and demand at nodes, respectively, and hence depend on the angular frequency,**

*Θ**ω*. The admittance matrix,

**, which contains the information regarding the topology of the system and its dynamics, determines the nodal outflows**

*Y***that are**

*Θ**admitted*from an input of nodal pressures

**. It can be evaluated (Zecchin**

*Ψ**et al.*2009) as where

*i*and

*k*denote two nodes of the system,

*λ*

_{j}is the

*j*-th branch of the system,

*Λ*

_{i}is the set of branches incident to node

*i*, and

*l*

_{j}is the length of the generic

*j*-th branch. In Equation (2), the propagation operator, , and the characteristic impedance, , can be evaluated for the

*j*-th link considering that

*α*

_{j}=

*ιω*[(

*ga*

_{j}/

*A*

^{2}

_{j}) + 2

*A*

_{j}

*f*

_{j}(

*ω*)] and

*β*

_{j}= [(

*ιω*/

*gA*

_{j}) + (

*f*

_{Fj}/2

*gD*

_{j}

*A*

^{2}

_{j})], where

*A*

_{j}is the pipe area,

*D*

_{j}is the diameter,

*f*

_{Fj}is the friction factor,

*a*

_{j}is the wave speed,

*g*is gravitational acceleration and . The term

*f*

_{j}(

*ω*) takes into account the pipe material constitutive law (Ferrante & Capponi 2017a, 2017b, 2017c). For an elastic pipe material

*f*

_{j}(

*ω*) = 0 while for a viscoelastic Standard Linear Solid model (SL), i.e. a series of a spring and a Kelvin–Voigt element, it is

*f*

_{j}(

*ω*) =

*k*

_{mj}(

*E*

_{Rj}+

*ιωη*

_{Rj})

^{−1}, where, for a thin-walled pipe of thickness

*e*

_{j},

*k*

_{mj}=

*f*

_{Rj}

*ρgD*

_{j}/(2

*e*

_{j}),

*E*

_{Rj}and

*η*

_{Rj}are the viscoelastic pipe material rheological parameters,

*ρ*is the water density and

*f*

_{Rj}is the restraint coefficient for the pipe (Chaudhry 2014).

Once the admittance matrix is inverted and the boundary conditions are given, the Fourier transform of the pressure head variation at the nodes can be simultaneously evaluated. The time domain corresponding variables, *H*(*t*) and *Q*(*t*), are then obtained by means of the inverse Fourier transform of ** Ψ**(

*ω*) and

**(**

*Θ**ω*), respectively.

The NAMM has been validated by means of numerical and laboratory tests on branched systems, considering viscoelastic and unsteady-friction effects, and used for diagnoses in terms of leak detection (Capponi *et al.* 2017a). Since a numerical investigation is presented, for the sake of simplicity the unsteady-friction effects are not implemented in the model.

## THE CONSIDERED SYSTEM

The numerical signals introduced in the following are generated by a NAMM model, considering the fast and complete closure of the valve V2 in the branched system of Figure 1(a). The system consists of three pipes, with internal diameter *D* = 93.3 mm and wall thickness *e* = 8.1 mm, of lengths *L*_{1} = 116.8 m, *L*_{2} = 61.8 m and *L*_{3} = 197.8 m, respectively. They connect the junction Y with the reservoir R, the valve V2 and the dead end V1, respectively. In the NAMM model, the node V2 is modeled as a demand-controlled node, i.e. the Fourier transform of the discharge variation is specified as a boundary condition and the pressure is considered as unknown.

The piezometric head time history at node V2 is simulated and then used as a benchmark signal in a calibration procedure where the lengths of the branches are considered as unknown. The calibration is based on the minimization of the mean squared error, *σ*^{2}, given by the sum of the squared differences between the benchmark and the simulated signals, divided by the number of the samples.

The same NAMM model is hence used twice, with different input and output. The direct solution of the transient equations, which yields the benchmark signal, uses the system geometry as input to determine the system response. Conversely, the calibration procedure aims at the solution of the inverse problem, considering the system response as the input and the branch lengths as unknowns.

Following a heuristic approach, in this preliminary numerical investigation all the possible combinations of the branch lengths of the system are considered. The configuration of the system, the boundary conditions and all the other geometric data are assumed as known. The first 20 s of the 100 Hz simulated signals are considered for the comparisons and hence 2,000 samples are used to evaluate the optimization function *σ*^{2}. Two pipe-material rheological behaviors are simulated, elastic and viscoelastic.

## CALIBRATION OF THE PIPE LENGTHS: ELASTIC MATERIAL

In the first simulation, the benchmark signal is generated at node V2 immediately upstream of the valve, assuming the same linear elastic material for all the pipes (Figure 2(a)). This assumption of the material rheology corresponds to using the constant wave speed *a* = 377.15 m/s for all the branches. The chosen value of the wave speed can be attributed, as an example, to oriented poly-vinyl chloride pipes, which exhibit a weak viscoelastic behavior, and it is much smaller than the values corresponding to the more common metallic pipes. This value is used here to enhance the comparison with the following viscoelastic pipe case, without any loss of significance for the shown results and conclusions.

### The shape of the optimization function

Since the number of the calibrated parameters is small and they can vary in a bounded and relatively small range, it is possible to take advantage of NAMM efficiency to analyze directly the variation of the optimization function *σ*^{2} on a grid in the parameter space. Considering that few seconds are required for the computation of each simulation, the computational burden for a grid of less than 10^{5} simulations is reasonable. This brute force calibration procedure, which allows the detection of the global minimum by the direct scrutiny of *σ*^{2} (Ferrante & Capponi 2017), is expected to be slower than more sophisticated techniques that need a reduced number of function evaluations and its accuracy depends on the grid discretization.

In a first set of calibrations, the link lengths are assumed unknown one at a time, while the other two are known. A grid for each branch length is defined between 0.0 m and 400.0 m with a step of Δ*l* = 3.0 m and *σ*^{2} is evaluated at each point of the grid. The discretization step is consistent with the length *a*/*f* ≈ 3.7 m corresponding to the travel distance of a pressure wave in the time period between two signal samples.

Direct scrutiny of the optimization function allows the detection of the global minimum for each branch length and, when the length value corresponds to the value used to generate the benchmark signal, it is *σ*^{2} = 0. The variations of *σ*^{2} with the lengths *L*_{1}, *L*_{2}, *L*_{3} are shown in Figure 3. It can be observed that, even if each figure clearly shows the global minimum of *σ*^{2}, there are other local minima in the analyzed range of variation. As an example, the local minimum for *L*_{3} = 35.4 m is particularly pronounced and many other minima are generated by the oscillation of *σ*^{2}. Due to the spatial grid step, the minimum for *L*_{2} corresponds to a value of *σ*^{2} greater than 10 m^{2}.

In the second set of calibration, two branch lengths are considered as unknown while the third is known. As a result, since two parameters are considered, *σ*^{2} is evaluated on a 2-D grid, as shown in Figure 4. The attained minima, denoted in Figure 4 by the white crosses, coincide with the true values on the branch lengths, with the given grid resolution. The function variation is irregular, with some V-shaped valleys containing local minima.

### Calibration by genetic and Nelder–Mead algorithms

To verify the effects of the optimization function shapes on the reliability of the calibration, a genetic algorithm (GA) and the Nelder–Mead NOA are used in series. The GA can explore a wider area of the parameter space, with respect to NOA, and locates several local minima. NOA explores a smaller area but can reach the global minimum with a better accuracy, with respect to GA. Hence the GA solution is used as the initial value for the NOA while the final results summarized in Table 1 are obtained by means of NOA.

No. . | Calibrated lengths . | Result [m] . | σ^{2} [m^{2}]
. | Success . |
---|---|---|---|---|

1 | L_{1} | 116.8 | 9.65·10^{−27} | yes |

2 | L_{2} | 61.8 | 9.06·10^{−27} | yes |

3 | L_{3} | 35.4 | 86.4 | no |

4 | L_{1};L_{2} | 116.8; 61.8 | 1.7·10^{−26} | yes |

5 | L_{1};L_{3} | 18.2; 17.4 | 86.2 | no |

6 | L_{2};L_{3} | 1.3; 2.1 | 83.3 | no |

7 | L_{1};L_{2};L_{3} | 2.0; 2.4; 0.1 | 85.4 | no |

No. . | Calibrated lengths . | Result [m] . | σ^{2} [m^{2}]
. | Success . |
---|---|---|---|---|

1 | L_{1} | 116.8 | 9.65·10^{−27} | yes |

2 | L_{2} | 61.8 | 9.06·10^{−27} | yes |

3 | L_{3} | 35.4 | 86.4 | no |

4 | L_{1};L_{2} | 116.8; 61.8 | 1.7·10^{−26} | yes |

5 | L_{1};L_{3} | 18.2; 17.4 | 86.2 | no |

6 | L_{2};L_{3} | 1.3; 2.1 | 83.3 | no |

7 | L_{1};L_{2};L_{3} | 2.0; 2.4; 0.1 | 85.4 | no |

GA requires an initial value of the parameter and a range defined by two bounds, a minimum and a maximum value. In the used optimization routine, an initial value for each parameter is assigned and then the routine evaluates the individual initial values by means of a random generation, considering the given bounds. The initial value for each length is set to 100 m for all the calibrations and the bounds for the GA are 0 m and 250 m. For each branch in the GA a population of 100 individuals is evolved through 50 generations, with an adaptive feasible mutation function and an intermediate crossover function, with a ratio of 0.5 and a crossover fraction of 0.8. For each parameter, the total number of 5,000 simulations evaluated by the GA is much greater than the number of solutions evaluated on the 135 nodes of the regular grid. The results of this procedure are summarized in Table 1.

In the last column a ‘yes’ distinguishes the successful calibrations that match the system geometry with the same given settings. Some calibrations are able to find the correct value of the length, even if there are two unknowns, such as in calibration no. 4. When the length *L*_{3} is one of the calibrated parameters, the calibration procedure fails (‘no’ in the last column of Table 1). The reasons for the failures can depend on the fact that the calibration procedure, with the given number of individuals, populations and generations, is not sophisticated enough to overcome problems that can be due to the shape of the optimization function. As an example, calibration no. 3 ends at the local minimum shown in Figure 3 for *L*_{3} = 35.4 m (the corresponding pressure signal is indicated by the dashed line in Figure 5).

The reasons for the failures could be due to the differences between the affected calibrations and the others, such as the fact that *L*_{3} is the longest length (and, furthermore, *L*_{3} > *L*_{1} + *L*_{2}) or that it is linked to the farthest node from the maneuver section where the signal is simulated. This means that, when the first wave reflected at the reservoir comes back to the maneuver section, there are a number of waves, reflected and transmitted at the other nodes of the system, that already have passed through the maneuver section. Due to the superposition of the effects of such waves it can be difficult to distinguish the wave reflected from the reservoir and to determine the length of the branch incident to it. The calibration failure could also depend on the type of boundary condition at node R, since the presence of a reservoir instead of other boundary conditions could affect the calibration of the length of the link incident to that node.

To investigate the causes of this failure and of the other failures of Table 1, a new system is considered, identical to the previous one except in the branches 2 and 3, to which the new lengths = *L*_{3} and = *L*_{2} are attributed, respectively (Figure 1(b)). The switch of the lengths makes it possible to investigate (1) if the calibration procedure can end successfully when the longest link is the closest to the maneuver section and (2) if it cannot find the length of the link incident to the reservoir node when its length is equal to *L*_{2}. A new benchmark signal is generated for the new system and the same procedure of the previous calibration is followed to determine and .

Using the same settings as the previous calibration, = 61.8 m is successfully found with *σ*^{2} = 1.9·10^{−27} m^{2}. This means that the boundary condition of the reservoir does not prejudice the success of the calibration. When = 197.8 m is calibrated, even if it is the branch incident to the maneuver node, it is not found by the chosen calibration procedure and the result is 0.03 m with *σ*^{2} = 119.9 m^{2}, which is an unreasonable solution.

The reasons for these results, summarized in Table 2, can be found in the shape of the optimization function for the lengths and shown in Figure 6.

No. . | Calibrated lengths . | Result [m] . | σ^{2} [m^{2}]
. | Success . |
---|---|---|---|---|

1 | 0.03 | 119.90 | no | |

2 | 61.8 | 1.90·10^{−27} | yes |

No. . | Calibrated lengths . | Result [m] . | σ^{2} [m^{2}]
. | Success . |
---|---|---|---|---|

1 | 0.03 | 119.90 | no | |

2 | 61.8 | 1.90·10^{−27} | yes |

It can be noted that the shape of *σ*^{2} in the case of the variation of is regular and does not present marked local minima that could mislead the calibration procedure. On the other hand, in the case of the variation of , the *σ*^{2} shape is jagged with several local minima and presents a sort of hump between = 0.0 m and = 150.0 m. Such a hump can be a significant obstacle for the calibration, leading the search for the solution to local minima and not to the global minimum.

This shape of the optimization function may explain the reduced reliability and the efficiency of the GA and NOA.

The algorithm settings are crucial for the success of the calibration because, for example, a proper choice of the initial value for the parameter to be calibrated makes the difference and allows the success of the calibration. In fact, in the case of , if the starting value is set to 100.0 m as in the previous calibrations, the algorithm finds a local minimum at one side of the hump. Otherwise, if the starting value is set, for example, to 150.0 m, the algorithm finds the global minimum with *σ*^{2} = 0. These observations indicate how important the knowledge is of the optimization function shape and the use of more sophisticated calibration, with larger computational effort.

## CALIBRATION OF THE PIPE LENGTHS: VISCOELASTIC MATERIAL

As in the elastic case, a numerical signal (Figure 2(b)) is generated for the system of Figure 1(a). The pipes are all made of the same viscoelastic material modeled by a Standard Linear Solid model with *a* = 377.15 m/s, *E*_{R} = 1.4·10^{10} Pa and *η*_{R} = 2.5·10^{9} Pa s.

### The shape of the optimization function

The shape of *σ*^{2} is firstly studied varying each one of the branch lengths at a time between 0.0 m and 400.0 m. The results are shown in Figure 7, which corresponds to Figure 3 for the elastic material case.

It can be noted that all these three curves are smoother than in the elastic case and this can be related to the damping caused by the viscoelasticity. Furthermore, when *L*_{3} is varied (Figure 7), *σ*^{2} is characterized by a shape with a hump that leads the optimization algorithm to the local minimum of 33.2 m of Table 3 (calibration no. 3), as in the elastic case. Also in this case, if the initial value for the optimization is set, as an example, to 150.0 m instead of 100.0 m, the calibration finds the exact value for *L*_{3}.

No. . | Calibrated parameters . | Result [m] . | σ^{2} [m^{2}]
. | Success . |
---|---|---|---|---|

1 | L_{1} | 116.8 | 1.1·10^{−27} | yes |

2 | L_{2} | 61.8 | 9.4·10^{−28} | yes |

3 | L_{3} | 33.2 | 35.7 | no |

4 | L_{1};L_{2} | 116.8; 61.8 | 7.9·10^{−28} | yes |

5 | L_{1};L_{3} | 119.2; 31.5 | 35.7 | no |

6 | L_{2};L_{3} | 61.8; 197.8 | 3.1·10^{−27} | yes |

7 | L_{1};L_{2};L_{3} | 173.0; 38.3; 87.7 | 33.2 | no |

No. . | Calibrated parameters . | Result [m] . | σ^{2} [m^{2}]
. | Success . |
---|---|---|---|---|

1 | L_{1} | 116.8 | 1.1·10^{−27} | yes |

2 | L_{2} | 61.8 | 9.4·10^{−28} | yes |

3 | L_{3} | 33.2 | 35.7 | no |

4 | L_{1};L_{2} | 116.8; 61.8 | 7.9·10^{−28} | yes |

5 | L_{1};L_{3} | 119.2; 31.5 | 35.7 | no |

6 | L_{2};L_{3} | 61.8; 197.8 | 3.1·10^{−27} | yes |

7 | L_{1};L_{2};L_{3} | 173.0; 38.3; 87.7 | 33.2 | no |

As a second step, the shape of *σ*^{2} is analyzed varying two of the branch lengths between 0.0 m and 400.0 m, considering the other as known. The results are shown in Figure 8, which corresponds to Figure 4 for the elastic material case. The global minima of the function are pointed out by white crosses and correspond to the true values on the branch lengths. As for the one-parameter analysis, the viscoelastic material produces a smoothing of *σ*^{2}.

### Calibration by genetic and Nelder–Mead algorithms

The generated benchmark signal is also used to estimate the system parameters using GA and NOA. The results of the calibrations are shown in Table 3, which corresponds to Table 1 for the elastic material pipes.

It can be noted that, as in the elastic material case, the algorithm can easily find the lengths of branches 1 and 2, while it cannot find the length of branch 3, *L*_{3}, alone or in combination with the other lengths. The only difference with respect to the elastic case is the combination of *L*_{3} with *L*_{2} (calibration no. 6), for which the NOA can find the correct values. This difference likely depends on the shape of the optimization function.

## CONCLUSIONS

The diagnosis of network systems by means of transients is generally based on the assumption that the lengths of the system branches are known. This requirement is not always fulfilled in practice, due to the lack of information about the actual system.

The models obtained by the integration of the transient governing equations in the frequency domain, or frequency domain models, explicitly use the length of the branches as parameters of the equations being solved. Hence, since usually the linearization needed for such an integration can be neglected or corrected, they can be implemented in a calibration procedure to estimate the length of the branches such as any other parameter.

In this paper the NAMM is used to generate a benchmark signal solving the transient direct problem in a branched system. Then, the signal is used as the input of an inverse problem to infer the length of the branches of the system by a calibration procedure.

The mean squared error is used in the calibration as the optimization function and the procedure is tested for both elastic and viscoelastic pipe materials, to assess the influence of pipe material rheology on the calibration results.

The shape of the optimization function is studied thanks to the model efficiency that allows the direct scrutiny of the solution in the parameter space, giving a heuristic interpretation of the existence of the solution. The function is investigated numerically on a regular grid, considering the variation of one parameter at a time or with two parameters simultaneously. In all the considered cases, the minima evaluated on the grid correspond to the true parameter values, i.e. to the branch lengths used for the benchmark signals, according to the used approximation. The differences between the estimated values and the true ones are always less than the used grid step, which is consistent with the time step and wave speed of the benchmark numerical signals.

The use of a GA to explore a wide range of the parameter variation and then of the Nelder–Mead optimization algorithm, to precisely point out the minimum, leads to the solution, but not for all the considered cases. The reasons can be explained by the used initial values and calibration settings and by the shape of the minimization function. The number of GA individuals, that is the number of the function evaluations for each generation, is fixed so that it is equal to the number of nodes of the grid where the function is also evaluated in the brute force calibration procedure.

In the viscoelastic case the results are slightly better than the elastic case, in terms of optimization algorithm, probably because the shape of the optimization function is smoother than the elastic case. This result indicates that viscoelasticity is not an obstacle in the diagnosis of networks, at least in terms of pipe length determination, but it facilitates the calibration in flattening the optimization function shape.

As a result of this work, it can be concluded that the NAMM can be applied to determine the length of the branches in a complex system. For the shown system, in the investigated range of the parameter variation, the solution exists and corresponds to a global minimum of the used optimization function. This result is not general but heuristic and limited to calibrations involving two pipe lengths as parameters, where the direct scrutiny of the function has been shown. This preliminary investigation can be considered as a first step toward more sophisticated optimization algorithms and more complex systems, where other nodes with unknown locations can be introduced. These nodes could also be leaks, changes in diameters, malfunctioning valves or other anomalies.

## ACKNOWLEDGMENTS

This research has been funded by the University of Perugia and by the Italian Ministry of Education, University and Research (MIUR) under the Project of Relevant National Interest ‘Tools and procedures for an advanced and sustainable management of water distribution systems’.

## REFERENCES

*(*American Water Works Association

*)*

*Laplace-Domain Analysis of Fluid Line Networks with Applications to Time-Domain Simulation and System Parameter Identification*. PhD Thesis, School of Civil, Environmental & Mining Engineering, Faculty of Engineering, Computer and Mathematical Sciences, the University of Adelaide, SA, Australia