- Email: [email protected]
Calibration of a model of an operational water distribution system containing pipes of different age
Advances in Engineering Software 40 (2009) 659–664
Contents lists available at ScienceDirect
Advances in Engineering Software journal homepage: www.elsevier.com/locate/advengsoft
Calibration of a model of an operational water distribution system containing pipes of different age T. Koppel, A. Vassiljev * Department of Mechanics, Tallinn University of Technology, Ehitajate tee 5, Tallinn 19086, Estonia
a r t i c l e
i n f o
Article history: Received 20 November 2007 Accepted 28 November 2008 Available online 20 February 2009 Keywords: Water distribution Calibration Roughness Levenberg–Marquardt Epanet
a b s t r a c t The aim of the paper is to demonstrate that the Levenberg–Marquardt algorithm can give successful results when operational water distribution systems are calibrated with the proper selection of parameter increment for the calculation of partial derivatives. The functional dependence of pipe roughness on age, which describes linear and nonlinear dependences, is proposed for the calibration of a model of a water distribution system containing pipes of different age. It is also shown that the visualization of response surface on a coarse grid is very useful for the analysis of the results. Ó 2009 Published by Elsevier Ltd.
1. Introduction All models of water distribution systems (WDS) require calibration. The precision of hydraulic models depends on how accurately they have been calibrated. The calibration parameters usually include pipe roughness, pipe diameters and demands, where the first two relate to flow conditions and the last one to boundary conditions. The corrosion and deposition processes, which occur over time after the pipe instalment, make it more difficult to determine the actual pipe diameter. Therefore, in the absence of another value, nominal pipe diameters are generally used for model development, and the roughness coefficient is adjusted to compensate the change in diameter due to the pipe wall build-up [1,2]. Economic developments in Estonia and the resulting increased water costs require that each consumer be supplied with a water flow meter, which in turn facilitates calibration. Calibration of pipe roughness is a difficult task if the WDS contains pipes of different age. In an ideal situation roughness coefficients for all pipes should be calibrated. WDS model calibration relies upon field measurement data (junction pressures, pipe flows, etc. [3]). However, measurements are quite expensive and usually much fewer measurements are available than required for calibration. Because of a large number of unknown values it is impossible to calibrate the model of a real system precisely [4]. Pipe roughness for all links cannot be determined accurately even when considerable efforts have been made to collect data [5]. Thus, a major problem associated with model calibration is the need to calibrate a * Corresponding author. E-mail address: [email protected] (A. Vassiljev). 0965-9978/$ - see front matter Ó 2009 Published by Elsevier Ltd. doi:10.1016/j.advengsoft.2008.11.015
large number of pipes using only a few measurements. Therefore, grouping of pipes is widely used at calibration [5]. All pipes are clustered into groups of identical estimated roughness and calibration is used to find these values. Roughness of old pipes depends on age, pipe material, water quality in WDS and some other factors. Clustering of pipes is a difficult task for old WDS because great age differences of pipes lead to a large number of clusters and consequently to a large number of parameters to be calibrated. Therefore, first the aim was to show that a flexible function, consisting of a small number of parameters to describe the dependence of roughness on age may be used instead of clustering. The parameters of the function may be found by an implicit method [6–8] that minimizes the objective function (OF), the sum of the squares of the differences between the measured and the simulated values. The simplest way to find the minimal value of the OF is the trialand-error algorithm. However, this algorithm would require a long period of time for all the possible variants to be examined. Several researchers have reported that standard algorithms (e.g. the Levenberg–Marquardt algorithm (LMA)) often find local minimums and they have proposed to use the Genetic Algorithm (GA) for calibration [9]. However, the GA also consumes much computer time. Kapelan et al. [10] proposed to use the combination of the GA and LMA algorithms in order to decrease computer time. A comprehensive review of available water distribution calibration models can be found in [11]. The calculations for an operational WDS containing thousands of pipes showed that the LMA had failed because of oscillations of the OF. The calculations of the OF indicated that its surface looks like emery paper with a large number of small local maximums and minimums. However, that surface would be very good for
660
T. Koppel, A. Vassiljev / Advances in Engineering Software 40 (2009) 659–664
the LMA (one minimum with good gradients) if we could eliminate these small variations. The simplest way to do it is to select a long enough step along the parameters, when calculating partial derivatives for the LMA [12]. This approach was tested and the second goal of the paper is to prove that the LMA algorithm works in this case quite well. The LMA requires only 1–2 min to find three parameters for a network containing more than 2000 pipes (computer with CPU 3.4 GHz). Short computer time will enable one to do more sophisticated analysis and to test many different variants of calibration. This investigation showed that the LMA is still useful for the calibration of operational WDS but the LMA needs information on the precision of calculations of the OF to work successfully. That is a significant result because the LMA requires much less computer time than the popular GA. The flexible function proposed in the article allows the dependence of roughness on age to be described for old WDS using a small number of parameters. It is shown also that the visualization of the OF surface is useful for the analysis of calibration results. 2. Method of calibration Optimization methods are used to search for a solution of how to describe the unknown calibration parameters which will minimize the OF [4,13], which usually is the sum of the residuals (the difference between the modelled and the observed values) squared. Since the number of measurements is usually much smaller than the number of pipes and nodes, it is clear that we cannot calibrate each single pipe. If a WDS contains pipes of different age, it will be reasonable to approximate the value of pipe roughness by the functional dependence on age and to try to find function parameters. Such approach will limit our problem to a search for a small number of the parameters of the function. Echávez [14] reports that the dependence of roughness on age is usually approximated by linear increment with time
e ¼ e0 þ at;
ð1Þ
where
e0 – roughness of the new pipe (mm), a – coefficient (mm/year), t – time (year).
Approximation (2) is quite flexible and describes both linear and nonlinear dependences with a different degree of nonlinearity. Eq. (1) is practically a particular case of Eq. (2) if b = 1.The same approximation for the Hazen–Williams C-factor will be
C i ¼ C min þ ðC max C min Þ½ðt max t i Þ=ðtmax tmin Þb ;
ð3Þ
where Ci – Hazen–Williams C-factor for pipe i. Thus, we have three parameters: emin, emax (or Cmin, Cmax) and the exponent b. If a WDS contains also new pipes, the number of parameters will be reduced to two because we know emin (or Cmax) for new pipes. The dependence of the roughness coefficient on a material may be approximated by assigning different emax for different materials. The use of Eqs. (2) and (3) is simpler than the clustering of pipes into several groups by age and finding different roughnesses for different groups if the WDS includes pipes of very different ages. Clustering results in too many groups in this case. Parameters emax, b (and emin if necessary) are found by minimizing the OF. Two methods were used for this purpose – the trialand-error method and the LMA algorithm. By the TOOLKIT developed by Rossman [15] for the EPANET2 calculations are made automatically. The TOOLKIT allows creation of external programs which manage the process of calculations and therefore it is very useful for calibration. Several subroutines were developed during this investigation using the TOOLKIT. One of them finds emax, emin and b for two groups of materials. Thus, three parameters for each of the two groups of the materials can be found. Each parameter may also be fixed if we know the value of the parameter in advance (for example, the roughness of brandnew pipes). Subroutines were developed in Visual Basic and MS Excel workbook was used for additional data. The TOOLKIT was also changed to work with double precision. The program includes software developed by the University of Chicago, as Operator of Argonne National Laboratory for nonlinear optimization, which was rewritten in Visual Basic and was modified to give the users a possibility to test the different minimal steps along the parameters for the calculation of partial derivatives in the LMA. The possibility to calculate standard errors of parameters was added to this software as well. 3. Results and discussion
The range of the coefficient a is 0.061–2.13 mm/year, depending on water quality, pipe material and other factors. Experiments [14] have also shown that the dependence of roughness on age may be nonlinear. Thus, experiments show that there is no universal dependence of roughness on age. Different water distribution systems have different dependence, which must be found practically for each model that describes the system with pipes of different age. Therefore, the type of functional dependence must be selected first. It would give us a possibility to find the linear function and the nonlinear function with a different degree of nonlinearity. For the pipe roughness e used to determine the friction factor for the Darcy– Weisbach formula it may be as follows:
ei ¼ emax ðemax emin Þ½ðtmax ti Þ=ðtmax tmin Þb ;
ð2Þ
where
ei – roughness for pipe i, mm, emin and emax – minimal and maximal roughness values which correspond to minimal and maximal ages of pipes, mm, tmin and tmax – minimal and maximal ages of pipes in the WDS, year.
A WDS representing a part of the whole WDS of the City of Tallinn (Fig. 1) was used as an example of the operational WDS. The characteristics of pipes installed are presented in Table 1. It can be seen that the age of pipes varies from 0 to 41 years. Cast iron is the material of the overwhelming majority of pipes (81% of the total length of pipes). Other materials used are plastic and steel (11% and 8% of the total length, respectively). Calibrations were performed on the basis of pressures measured during nine days (from wednesday of the first week to thursday of the following week). Pressures were measured in 25 uniformly distributed points (Fig. 1). The time interval of measurements was 15 min. Water flow was measured in two points of water inlets into the system with the time interval of 1 min (Fig. 1). The WDS includes over 1000 consumers, each supplied with a flow meter. The consumption is read once a month. These data enable one to calculate the part of the total demand consumed by each user. As the consumption is not constant during the day, the diurnal dynamics was investigated for different groups of consumers. All consumers were divided into four groups: multi-storey apartment buildings; small apartment buildings; large social buildings; small social buildings. The diurnal dynamics of water
661
T. Koppel, A. Vassiljev / Advances in Engineering Software 40 (2009) 659–664
Fig. 1. Modelled Network Geometry (black dots indicate the points of pressure measurement, black squares are the points of water inflow (pressure and flow rate measurements)).
Table 1 Age, material, total length and diameter of pipes in the WDS. Ages
Pipe material Cast iron
Plastic
Steel
Total length, m
Mean diameter, mm
Total length, m
Mean diameter, mm
0–1 2–3 4–6 19–21 24–25 26–27 28–30 31–32 33–37 38–41
10 490 2048 470 13,556 6740 1783 35,223 4370 5580
150 230 166 150 170 200 335 206 185 140
2043 3000 4245
160 155 109
Total
70,270
9287
demand by consumers of different groups was measured for several consumers of each group. Then the dynamics was applied to all other consumers of the same group. Leakages were estimated using the minimal night flow method [16]. According to flow measurements, 82% of water is consumed by multi-storey apartment buildings, 2% by small buildings, 4% by social buildings (e.g. schools, hospitals), and 1% by other small social buildings and 8% by leakages. We assume that the remaining 3% (or unaccounted demand) was used by multi-storey apartment buildings. The assumption was made on the basis of the analysis of diurnal dynamics of water inflow into the system, the dynamics of consumption by different consumers and the differences between water inflow and total consumption. The dynamics of differences was very close to the dynamics of consumption by multi-storey apartment buildings. Therefore, the differences between water inflow and demand were distributed between multi-storey apartment buildings proportional to the demand. Leakages were distributed using the approach proposed in [17].
Total length, m
Mean diameter, mm
2832 375 1150 2350
458 360 330 500
7317
Low velocities in the WDS pipes of the City of Tallinn made the calibration process complicated. Meier and Barkdoll [18] recommend working with velocities 0.3 m/s or higher. The analysis of velocities showed that the network meets the above condition only at the highest demand. Walski [4] recommends calibration to be performed at maximal head loss (practically at the highest demand in our case). As a result of data analysis, two days with the highest demand and containing data on all loggers without missing data periods were selected from the set of data. One day was selected for calibration and the other for validation. The calibration process attempts to adjust model parameters in such a way that the field observations and the predicted results will be in a reasonable agreement. The objective function (4) was used to evaluate the agreement
OF ¼
n X i¼1
ðHoi Hpi Þ2 ;
ð4Þ
662
T. Koppel, A. Vassiljev / Advances in Engineering Software 40 (2009) 659–664
where OF – objective function to be minimized, Hoi and Hpi – observed and predicted pressure in the node i, n – number of pressure observations (25 in our case). All pipes were clustered into two groups: metal pipes and plastic pipes. Eq. (2) was used to present the dependence of roughness on age. Since the system contains new pipes as well (pipes of age = 0), it was assumed that emin is equal to the roughness of new pipes (0.26 mm for metal pipes, and 0.002 for plastic pipes). Our first attempt to find emax and exponent b by the LMA algorithm was unsuccessful. The results were always very close to the initial values of the parameters. After that we used the trial-and-error algorithm to investigate the response surface of the OF. Exponent b was assumed to equal 1 for both groups of pipes to simplify calculations. Thus, it was necessary to find two parameters: emax for metal pipes and emax for plastic pipes. The parameter domain was split into the coarse grid of values and the OF was calculated for each set of parameters. The obtained data were used to construct the response surface (Fig. 2). The minimal value of the OF is at the point where emax is equal to 14 mm for metal pipes and 0.002 mm for plastic pipes. It means that the maximal roughness (pipe roughness of the age of 41 years) is 14 mm for metal pipes. The roughness of plastic pipes had not changed within 6 years. The attempt to use the LMA algorithm was unsuccessful for this simple case as well. It was surprising because the response surface (Fig. 2) is a long narrow valley and the initial point was selected on the side of the valley. However, the LMA algorithm gave us again values which were very close to the initial values. Therefore, the calculations of the OF were made on the fine grid around the point on the side of the valley. The point with the coordinates of 15 mm for metal pipes and 0.2 mm for plastic pipes was selected as the central point. Fig. 2 shows that this point is located on the obvious slope of the response surface. The step of fine grid was selected to be 0.0001 mm for both axes. The results of calculations are presented in Fig. 3. It can be seen that the surface has a large number of local minimums. It is most likely that such behaviour is the result of the iterative technique used in the Epanet. The reason of the phenomenon may not be clear, but it does exist. The calculations also showed that the OF surface is often not smooth even for a simple WDS (10–20 pipes) if we use the Darcy–Weisbach formula to compute a head loss.
Fig. 2. Dependence of the objective function on the maximal roughness of metal and plastic pipes.
Fig. 3. Objective function on the fine grid.
The surface for a simple WDS is quite smooth if we use the Hazen–Williams formula, but it was not smooth for the network under consideration. Thus, the analysis showed that one of the reasons, that prevents us from using the LMA algorithm for calibration of the operational WDS is the non-smooth response surface of the OF. It means that the number of reliable digits returned by the OF must be taken into account [11]. This number is used to calculate the finite difference step size. Fig. 2 shows that the surface is not complex and the LMA algorithm can be used, but Fig. 3 shows that the finite difference step size must be sufficiently great for the LMA algorithm not to be sensitive to the oscillations of the OF. This step is estimated on the basis of the number of reliable digits returned by the OF. The number of reliable digits may be estimated by Hamming’s method [19] described in [20]. This number estimated by Hamming’s method was used to calculate the finite difference step size, which was successfully used in the LMA method. However, the procedure for the estimation of reliable digits and the step is not simple and needs additional calculations. A simpler way is to estimate the size of the step by trying out different step sizes for calculations. The step size appropriate for the LMA algorithm that gives the same results with different initial values of parameters was assumed as the usable step size [1]. The possibility to test different sizes was added to the software [21]. In our case, the LMA algorithm worked when the step was larger than 0.01% of the parameter value. A similar approach has also been used in [22]. The modified LMA algorithm showed also (as the trial-and-error algorithm) that the roughness of plastic pipes did not change. As a result, an attempt was made to find emax and exponents for the groups of cast iron and steel pipes (assuming that the roughness of plastic pipes remains unchanged and is equal to the roughness of new pipes 0.002 mm). The calculations showed that emax for cast iron is five times larger than for steel pipes, but the standard errors of emax and exponents were very high (several times higher than the parameters themselves). The next attempt was to estimate emax for these two groups, assuming that the exponents are equal to 1. However, the standard error for the maximal roughness of steel pipes was very great even for this simple case. High errors mean that the number of parameters is too large for our data set. Therefore, the initial variant of pipe grouping was used: metal pipes and plastic pipes. We assumed the roughness for plastic pipes to be 0.002 mm. An attempt was made to find emax and exponent b for the group of metal pipes. The results were as follows: emax = 15.317 mm (standard error = 0.973), b = 0.499 (standard error = 0.178). It can be seen that the error for exponent b was still
663
T. Koppel, A. Vassiljev / Advances in Engineering Software 40 (2009) 659–664
18 16 Roughness, mm
quite high. Therefore, the trial-and-error method was again used to construct the surface of the OF. The parameters obtained by the LMA (emax = 15.317, b = 0.499) were selected as the central point. The other points were calculated by the increment and decrement of the central point values by a fixed step. The step for emax was assumed to be 0.5 mm and for exponent 0.05. The number of steps was 10 for each parameter decrement and 9 for each parameter increment. Fig. 4 shows the response surface for the grid 20 20, representing 400 calculations by the Epanet2. It can be seen that the response surface is quite smooth with a thalweg – a line, connecting points of the surface with low values of the OF in the valley. The OF values are almost the same along the thalweg. Thus, practically it is possible to find exponent b for each emax, which gives almost the same minimal OF. The higher the emax, the lower the exponent b is. For example, the OF is equal to 6.80 for emax = 15.317 and b = 0.499 (central point) and 6.81 for emax = 10.3 and b = 1.0. Fig. 5 presents two curves – one with emax = 15.3 (corresponding to b = 0.5) and the other one with lower emax = 10.3 (corresponding to b = 1.0). The differences between nonlinear and linear variants are not great for age in the range of 0–30 years and practically coincide for the age of 30–32 years. It is not surprising because the overwhelming majority of pipes are of the age of 31–32 years (Table 1). However, the great differences for the age of 38–41 years are quite surprising, because the total length of pipes with this age is considerable (Table 1). The analysis of the location of pipes with different age, points of measurement and velocities showed that the pressure was measured in one point only in the oldest part of the WDS and the velocities in the area were very low (0.05 m/ s) even at the highest consumption. That is why both curves (1 and 2 in Fig. 5) give practically the same values of the OF and the surface of the OF has such a shape with thalweg (see Fig. 4). It means that the set of measurements is not sufficient for the evaluation of the exponent with a fair degree of confidence. The points of measurement must be selected by more complicated techniques than in the case of uniform distribution. Therefore, a decision was made to use the linear case of Eq. (2) (exponent b was assumed to be equal to 1.0), to find emax for this variant and use the variant as
1)
14
2)
12 10 8 6 4 2 0
0
10
20 Age, years
30
40
Fig. 5. Dependence of pipe roughness on the age of the pipe with different sets of parameters: (1) emax = 15.3, b = 0.5, OF = 6.80 and (2) emax = 10.3, b = 1.0, OF = 6.81.
the result of calibration. The resulting emax obtained by the LMA with the fixed exponent value b = 1.0 is equal to 10.4 with a standard error of 0.49. Thus, the analysis showed that the visualization of response surface is very useful to understand the results obtained by the LMA. The preparation of such a surface does not take much computer time on the coarse grid (e.g. a grid containing 21 values for each parameter needs 441 calculations by the Epanet2 and takes 4 min on the computer with CPU 3.4 GHz). As was mentioned above, the measurements were divided into two parts. The first part was used to calibrate the model. The second part was used for validation, i.e. the parameters obtained by calibration in the first part of series were used in the calculations for the second part of series, to verify model performance on the independent data. The calculations showed that the sum of squares of differences between the measured and the calculated values was almost the same for the validation of data (the value of the OF was 6.82 for the validation data, as compared to the value 6.80 for the part of data used for calibration). The comparison of the observed and the simulated pressures (validation data) is presented in Fig. 6. The points corresponding to the minimal demand are lower than the other points because at night the pressure decreased in the system.
70
Pressure simulated, m
60
Mean demand Maximal demand
50
Minimal demand
40 30 20 10 0 0
10
20
30
40
50
60
70
Pressure observed, m
Fig. 4. Surface of the objective function on the grid 20 20 points.
Fig. 6. Simulated pressure vs. observed pressure for mean, maximal and minimal demands.
664
T. Koppel, A. Vassiljev / Advances in Engineering Software 40 (2009) 659–664
4. Conclusions The Levenberg–Marquardt algorithm requires a low limit for the step along parameters to calculate partial derivatives in order to work on the operational WDS model containing thousands of pipes. The approximation of the dependence of pipe roughness on age by the function with a low number of parameters is easier than pipe grouping according to age. The analysis showed that the response surface is very useful to understand the results. The combination of the Levenberg–Marquardt algorithm with the response surface obtained by the trial-and-error method facilitates calibration significantly.
Acknowledgements This work was financially supported by Target Financing of Estonia (grant SF0140072s08) at Tallinn University of Technology and by the Estonian Science Foundation (grant G7646). References [1] Vassiljev A, Koppel T. Calibration of the model of an operational water distribution system. In: Topping BHV, Montero G, Montenegro R, editors. Proceedings of the fifth international conference on engineering computational technology. Stirlingshire, United Kingdom: Civil-Comp Press; 2006 [paper 155]. [2] Lansey KE, El-Shorbagy W, Ahmed I, Araujo J, Haan CT. Calibration assessment and data collection for water distribution networks. J Hydraul Eng 2001;127:270–9. [3] Greco M, Giudice GD. New approach to water distribution network calibration. J Hydraul Eng 1999;125(8):849–54. [4] Walski TM. Model calibration data: the good, the bad, and the useless. J AWWA 2000;92(1):94–9. [5] Mallick KN, Ahmed I, Tickle KS, Lansey KE. Determining pipe groupings for water distribution networks. J Water Resour Plan Manage, ASCE 2002;128(2):130–9.
[6] Ormsbee LE. Implicit network calibration. J Water Resour Plan Manage, ASCE 1989;115(2):243–57. [7] Lansey KE, Basnet C. Parameter estimation for water distribution networks. J Water Resour Plan Manage, ASCE 1991;117(1):126–44. [8] Liggett JA, Chen LC. Inverse transient analysis in pipe networks. J Hydraul Eng, ASCE 1994;120(8):934. [9] Vitkovsky JP, Simpson AR, Lambert MF. Leak detection and calibration using transients and genetic algorithms. J Water Resour Plan Manage, ASCE 2000;126(4):262. [10] Kapelan ZS, Savic DA, Walters GA. A hybrid inverse transient model for leakage detection and roughness calibration in pipe networks. J Hydraul Res 2003;41(5):481–92. [11] Methods Haestad, Walski TM, Chase DV, Savic DA, Grayman WM, Beckwith S, et al. Advanced water distribution modeling and management. Waterbury, USA: Haestead Press; 2003. [12] Dennis Jr JE, Schnabel Robert B. Numerical methods for unconstrained optimization and nonlinear equations. Philadelphia: SIAM; 1996. [13] Reddy PVN, Sridharan K, Rao PV. WLS method for parameter estimation in water distribution networks. J Water Resour Plan Manage, ASCE 1996;122(3):157. [14] Echávez G. Increase in losses coefficient with age for small diameter pipes. J Hydraul Eng 1997;2:157–9. [15] Rossman LA. EPANET 2 User’s Manual. Water Supply and Water Resources Division, National Risk Management Research Laboratory, Cincinnati, OH 45268, September 2002. [16] Garcia VJ, Cabrera E, Cabrera E, Jr. The minimum night flow method revisited. In: 8th annual water distribution systems analysis symposium, Cincinnati, Ohio, USA, 27–30 August 2006. [17] Koppel T, Kändler T, Vassiljev A. Calibration of the water distribution network model. In: Brebbia CA, Anagnostopoulos P, Katsifarakis K, Cheng H-D, editors. Proceedings of the water resources management conference, Halkidiki, Greece, September 24–26, 2001. WIT Press: Southampton, Boston; 2001. p. 207–17. [18] Meier RW, Barkdoll BD. Sampling. Design for network model calibration using genetic algorithms. J Water Resour Plan Manage, ASCE 2000;126(4):245–50. [19] Hamming RW. Numerical methods for scientists. 2nd ed. New York: McGrawHill Book Co.; 1973. [20] Gill E, Murray W, Wright M. Practical optimization. New York: Academic Press; 1981. [21] Koppel T, Vassiljev A. Experience in modelling of water supply systems. Congress Publications of ECWATECH-2004, 1–4 June, 2004. [22] Deepan B. Application of optimization methods to calibration of water distribution systems. In: World water and environmental resources congress, June 27–July 1, 2004.
Contents lists available at ScienceDirect
Advances in Engineering Software journal homepage: www.elsevier.com/locate/advengsoft
Calibration of a model of an operational water distribution system containing pipes of different age T. Koppel, A. Vassiljev * Department of Mechanics, Tallinn University of Technology, Ehitajate tee 5, Tallinn 19086, Estonia
a r t i c l e
i n f o
Article history: Received 20 November 2007 Accepted 28 November 2008 Available online 20 February 2009 Keywords: Water distribution Calibration Roughness Levenberg–Marquardt Epanet
a b s t r a c t The aim of the paper is to demonstrate that the Levenberg–Marquardt algorithm can give successful results when operational water distribution systems are calibrated with the proper selection of parameter increment for the calculation of partial derivatives. The functional dependence of pipe roughness on age, which describes linear and nonlinear dependences, is proposed for the calibration of a model of a water distribution system containing pipes of different age. It is also shown that the visualization of response surface on a coarse grid is very useful for the analysis of the results. Ó 2009 Published by Elsevier Ltd.
1. Introduction All models of water distribution systems (WDS) require calibration. The precision of hydraulic models depends on how accurately they have been calibrated. The calibration parameters usually include pipe roughness, pipe diameters and demands, where the first two relate to flow conditions and the last one to boundary conditions. The corrosion and deposition processes, which occur over time after the pipe instalment, make it more difficult to determine the actual pipe diameter. Therefore, in the absence of another value, nominal pipe diameters are generally used for model development, and the roughness coefficient is adjusted to compensate the change in diameter due to the pipe wall build-up [1,2]. Economic developments in Estonia and the resulting increased water costs require that each consumer be supplied with a water flow meter, which in turn facilitates calibration. Calibration of pipe roughness is a difficult task if the WDS contains pipes of different age. In an ideal situation roughness coefficients for all pipes should be calibrated. WDS model calibration relies upon field measurement data (junction pressures, pipe flows, etc. [3]). However, measurements are quite expensive and usually much fewer measurements are available than required for calibration. Because of a large number of unknown values it is impossible to calibrate the model of a real system precisely [4]. Pipe roughness for all links cannot be determined accurately even when considerable efforts have been made to collect data [5]. Thus, a major problem associated with model calibration is the need to calibrate a * Corresponding author. E-mail address: [email protected] (A. Vassiljev). 0965-9978/$ - see front matter Ó 2009 Published by Elsevier Ltd. doi:10.1016/j.advengsoft.2008.11.015
large number of pipes using only a few measurements. Therefore, grouping of pipes is widely used at calibration [5]. All pipes are clustered into groups of identical estimated roughness and calibration is used to find these values. Roughness of old pipes depends on age, pipe material, water quality in WDS and some other factors. Clustering of pipes is a difficult task for old WDS because great age differences of pipes lead to a large number of clusters and consequently to a large number of parameters to be calibrated. Therefore, first the aim was to show that a flexible function, consisting of a small number of parameters to describe the dependence of roughness on age may be used instead of clustering. The parameters of the function may be found by an implicit method [6–8] that minimizes the objective function (OF), the sum of the squares of the differences between the measured and the simulated values. The simplest way to find the minimal value of the OF is the trialand-error algorithm. However, this algorithm would require a long period of time for all the possible variants to be examined. Several researchers have reported that standard algorithms (e.g. the Levenberg–Marquardt algorithm (LMA)) often find local minimums and they have proposed to use the Genetic Algorithm (GA) for calibration [9]. However, the GA also consumes much computer time. Kapelan et al. [10] proposed to use the combination of the GA and LMA algorithms in order to decrease computer time. A comprehensive review of available water distribution calibration models can be found in [11]. The calculations for an operational WDS containing thousands of pipes showed that the LMA had failed because of oscillations of the OF. The calculations of the OF indicated that its surface looks like emery paper with a large number of small local maximums and minimums. However, that surface would be very good for
660
T. Koppel, A. Vassiljev / Advances in Engineering Software 40 (2009) 659–664
the LMA (one minimum with good gradients) if we could eliminate these small variations. The simplest way to do it is to select a long enough step along the parameters, when calculating partial derivatives for the LMA [12]. This approach was tested and the second goal of the paper is to prove that the LMA algorithm works in this case quite well. The LMA requires only 1–2 min to find three parameters for a network containing more than 2000 pipes (computer with CPU 3.4 GHz). Short computer time will enable one to do more sophisticated analysis and to test many different variants of calibration. This investigation showed that the LMA is still useful for the calibration of operational WDS but the LMA needs information on the precision of calculations of the OF to work successfully. That is a significant result because the LMA requires much less computer time than the popular GA. The flexible function proposed in the article allows the dependence of roughness on age to be described for old WDS using a small number of parameters. It is shown also that the visualization of the OF surface is useful for the analysis of calibration results. 2. Method of calibration Optimization methods are used to search for a solution of how to describe the unknown calibration parameters which will minimize the OF [4,13], which usually is the sum of the residuals (the difference between the modelled and the observed values) squared. Since the number of measurements is usually much smaller than the number of pipes and nodes, it is clear that we cannot calibrate each single pipe. If a WDS contains pipes of different age, it will be reasonable to approximate the value of pipe roughness by the functional dependence on age and to try to find function parameters. Such approach will limit our problem to a search for a small number of the parameters of the function. Echávez [14] reports that the dependence of roughness on age is usually approximated by linear increment with time
e ¼ e0 þ at;
ð1Þ
where
e0 – roughness of the new pipe (mm), a – coefficient (mm/year), t – time (year).
Approximation (2) is quite flexible and describes both linear and nonlinear dependences with a different degree of nonlinearity. Eq. (1) is practically a particular case of Eq. (2) if b = 1.The same approximation for the Hazen–Williams C-factor will be
C i ¼ C min þ ðC max C min Þ½ðt max t i Þ=ðtmax tmin Þb ;
ð3Þ
where Ci – Hazen–Williams C-factor for pipe i. Thus, we have three parameters: emin, emax (or Cmin, Cmax) and the exponent b. If a WDS contains also new pipes, the number of parameters will be reduced to two because we know emin (or Cmax) for new pipes. The dependence of the roughness coefficient on a material may be approximated by assigning different emax for different materials. The use of Eqs. (2) and (3) is simpler than the clustering of pipes into several groups by age and finding different roughnesses for different groups if the WDS includes pipes of very different ages. Clustering results in too many groups in this case. Parameters emax, b (and emin if necessary) are found by minimizing the OF. Two methods were used for this purpose – the trialand-error method and the LMA algorithm. By the TOOLKIT developed by Rossman [15] for the EPANET2 calculations are made automatically. The TOOLKIT allows creation of external programs which manage the process of calculations and therefore it is very useful for calibration. Several subroutines were developed during this investigation using the TOOLKIT. One of them finds emax, emin and b for two groups of materials. Thus, three parameters for each of the two groups of the materials can be found. Each parameter may also be fixed if we know the value of the parameter in advance (for example, the roughness of brandnew pipes). Subroutines were developed in Visual Basic and MS Excel workbook was used for additional data. The TOOLKIT was also changed to work with double precision. The program includes software developed by the University of Chicago, as Operator of Argonne National Laboratory for nonlinear optimization, which was rewritten in Visual Basic and was modified to give the users a possibility to test the different minimal steps along the parameters for the calculation of partial derivatives in the LMA. The possibility to calculate standard errors of parameters was added to this software as well. 3. Results and discussion
The range of the coefficient a is 0.061–2.13 mm/year, depending on water quality, pipe material and other factors. Experiments [14] have also shown that the dependence of roughness on age may be nonlinear. Thus, experiments show that there is no universal dependence of roughness on age. Different water distribution systems have different dependence, which must be found practically for each model that describes the system with pipes of different age. Therefore, the type of functional dependence must be selected first. It would give us a possibility to find the linear function and the nonlinear function with a different degree of nonlinearity. For the pipe roughness e used to determine the friction factor for the Darcy– Weisbach formula it may be as follows:
ei ¼ emax ðemax emin Þ½ðtmax ti Þ=ðtmax tmin Þb ;
ð2Þ
where
ei – roughness for pipe i, mm, emin and emax – minimal and maximal roughness values which correspond to minimal and maximal ages of pipes, mm, tmin and tmax – minimal and maximal ages of pipes in the WDS, year.
A WDS representing a part of the whole WDS of the City of Tallinn (Fig. 1) was used as an example of the operational WDS. The characteristics of pipes installed are presented in Table 1. It can be seen that the age of pipes varies from 0 to 41 years. Cast iron is the material of the overwhelming majority of pipes (81% of the total length of pipes). Other materials used are plastic and steel (11% and 8% of the total length, respectively). Calibrations were performed on the basis of pressures measured during nine days (from wednesday of the first week to thursday of the following week). Pressures were measured in 25 uniformly distributed points (Fig. 1). The time interval of measurements was 15 min. Water flow was measured in two points of water inlets into the system with the time interval of 1 min (Fig. 1). The WDS includes over 1000 consumers, each supplied with a flow meter. The consumption is read once a month. These data enable one to calculate the part of the total demand consumed by each user. As the consumption is not constant during the day, the diurnal dynamics was investigated for different groups of consumers. All consumers were divided into four groups: multi-storey apartment buildings; small apartment buildings; large social buildings; small social buildings. The diurnal dynamics of water
661
T. Koppel, A. Vassiljev / Advances in Engineering Software 40 (2009) 659–664
Fig. 1. Modelled Network Geometry (black dots indicate the points of pressure measurement, black squares are the points of water inflow (pressure and flow rate measurements)).
Table 1 Age, material, total length and diameter of pipes in the WDS. Ages
Pipe material Cast iron
Plastic
Steel
Total length, m
Mean diameter, mm
Total length, m
Mean diameter, mm
0–1 2–3 4–6 19–21 24–25 26–27 28–30 31–32 33–37 38–41
10 490 2048 470 13,556 6740 1783 35,223 4370 5580
150 230 166 150 170 200 335 206 185 140
2043 3000 4245
160 155 109
Total
70,270
9287
demand by consumers of different groups was measured for several consumers of each group. Then the dynamics was applied to all other consumers of the same group. Leakages were estimated using the minimal night flow method [16]. According to flow measurements, 82% of water is consumed by multi-storey apartment buildings, 2% by small buildings, 4% by social buildings (e.g. schools, hospitals), and 1% by other small social buildings and 8% by leakages. We assume that the remaining 3% (or unaccounted demand) was used by multi-storey apartment buildings. The assumption was made on the basis of the analysis of diurnal dynamics of water inflow into the system, the dynamics of consumption by different consumers and the differences between water inflow and total consumption. The dynamics of differences was very close to the dynamics of consumption by multi-storey apartment buildings. Therefore, the differences between water inflow and demand were distributed between multi-storey apartment buildings proportional to the demand. Leakages were distributed using the approach proposed in [17].
Total length, m
Mean diameter, mm
2832 375 1150 2350
458 360 330 500
7317
Low velocities in the WDS pipes of the City of Tallinn made the calibration process complicated. Meier and Barkdoll [18] recommend working with velocities 0.3 m/s or higher. The analysis of velocities showed that the network meets the above condition only at the highest demand. Walski [4] recommends calibration to be performed at maximal head loss (practically at the highest demand in our case). As a result of data analysis, two days with the highest demand and containing data on all loggers without missing data periods were selected from the set of data. One day was selected for calibration and the other for validation. The calibration process attempts to adjust model parameters in such a way that the field observations and the predicted results will be in a reasonable agreement. The objective function (4) was used to evaluate the agreement
OF ¼
n X i¼1
ðHoi Hpi Þ2 ;
ð4Þ
662
T. Koppel, A. Vassiljev / Advances in Engineering Software 40 (2009) 659–664
where OF – objective function to be minimized, Hoi and Hpi – observed and predicted pressure in the node i, n – number of pressure observations (25 in our case). All pipes were clustered into two groups: metal pipes and plastic pipes. Eq. (2) was used to present the dependence of roughness on age. Since the system contains new pipes as well (pipes of age = 0), it was assumed that emin is equal to the roughness of new pipes (0.26 mm for metal pipes, and 0.002 for plastic pipes). Our first attempt to find emax and exponent b by the LMA algorithm was unsuccessful. The results were always very close to the initial values of the parameters. After that we used the trial-and-error algorithm to investigate the response surface of the OF. Exponent b was assumed to equal 1 for both groups of pipes to simplify calculations. Thus, it was necessary to find two parameters: emax for metal pipes and emax for plastic pipes. The parameter domain was split into the coarse grid of values and the OF was calculated for each set of parameters. The obtained data were used to construct the response surface (Fig. 2). The minimal value of the OF is at the point where emax is equal to 14 mm for metal pipes and 0.002 mm for plastic pipes. It means that the maximal roughness (pipe roughness of the age of 41 years) is 14 mm for metal pipes. The roughness of plastic pipes had not changed within 6 years. The attempt to use the LMA algorithm was unsuccessful for this simple case as well. It was surprising because the response surface (Fig. 2) is a long narrow valley and the initial point was selected on the side of the valley. However, the LMA algorithm gave us again values which were very close to the initial values. Therefore, the calculations of the OF were made on the fine grid around the point on the side of the valley. The point with the coordinates of 15 mm for metal pipes and 0.2 mm for plastic pipes was selected as the central point. Fig. 2 shows that this point is located on the obvious slope of the response surface. The step of fine grid was selected to be 0.0001 mm for both axes. The results of calculations are presented in Fig. 3. It can be seen that the surface has a large number of local minimums. It is most likely that such behaviour is the result of the iterative technique used in the Epanet. The reason of the phenomenon may not be clear, but it does exist. The calculations also showed that the OF surface is often not smooth even for a simple WDS (10–20 pipes) if we use the Darcy–Weisbach formula to compute a head loss.
Fig. 2. Dependence of the objective function on the maximal roughness of metal and plastic pipes.
Fig. 3. Objective function on the fine grid.
The surface for a simple WDS is quite smooth if we use the Hazen–Williams formula, but it was not smooth for the network under consideration. Thus, the analysis showed that one of the reasons, that prevents us from using the LMA algorithm for calibration of the operational WDS is the non-smooth response surface of the OF. It means that the number of reliable digits returned by the OF must be taken into account [11]. This number is used to calculate the finite difference step size. Fig. 2 shows that the surface is not complex and the LMA algorithm can be used, but Fig. 3 shows that the finite difference step size must be sufficiently great for the LMA algorithm not to be sensitive to the oscillations of the OF. This step is estimated on the basis of the number of reliable digits returned by the OF. The number of reliable digits may be estimated by Hamming’s method [19] described in [20]. This number estimated by Hamming’s method was used to calculate the finite difference step size, which was successfully used in the LMA method. However, the procedure for the estimation of reliable digits and the step is not simple and needs additional calculations. A simpler way is to estimate the size of the step by trying out different step sizes for calculations. The step size appropriate for the LMA algorithm that gives the same results with different initial values of parameters was assumed as the usable step size [1]. The possibility to test different sizes was added to the software [21]. In our case, the LMA algorithm worked when the step was larger than 0.01% of the parameter value. A similar approach has also been used in [22]. The modified LMA algorithm showed also (as the trial-and-error algorithm) that the roughness of plastic pipes did not change. As a result, an attempt was made to find emax and exponents for the groups of cast iron and steel pipes (assuming that the roughness of plastic pipes remains unchanged and is equal to the roughness of new pipes 0.002 mm). The calculations showed that emax for cast iron is five times larger than for steel pipes, but the standard errors of emax and exponents were very high (several times higher than the parameters themselves). The next attempt was to estimate emax for these two groups, assuming that the exponents are equal to 1. However, the standard error for the maximal roughness of steel pipes was very great even for this simple case. High errors mean that the number of parameters is too large for our data set. Therefore, the initial variant of pipe grouping was used: metal pipes and plastic pipes. We assumed the roughness for plastic pipes to be 0.002 mm. An attempt was made to find emax and exponent b for the group of metal pipes. The results were as follows: emax = 15.317 mm (standard error = 0.973), b = 0.499 (standard error = 0.178). It can be seen that the error for exponent b was still
663
T. Koppel, A. Vassiljev / Advances in Engineering Software 40 (2009) 659–664
18 16 Roughness, mm
quite high. Therefore, the trial-and-error method was again used to construct the surface of the OF. The parameters obtained by the LMA (emax = 15.317, b = 0.499) were selected as the central point. The other points were calculated by the increment and decrement of the central point values by a fixed step. The step for emax was assumed to be 0.5 mm and for exponent 0.05. The number of steps was 10 for each parameter decrement and 9 for each parameter increment. Fig. 4 shows the response surface for the grid 20 20, representing 400 calculations by the Epanet2. It can be seen that the response surface is quite smooth with a thalweg – a line, connecting points of the surface with low values of the OF in the valley. The OF values are almost the same along the thalweg. Thus, practically it is possible to find exponent b for each emax, which gives almost the same minimal OF. The higher the emax, the lower the exponent b is. For example, the OF is equal to 6.80 for emax = 15.317 and b = 0.499 (central point) and 6.81 for emax = 10.3 and b = 1.0. Fig. 5 presents two curves – one with emax = 15.3 (corresponding to b = 0.5) and the other one with lower emax = 10.3 (corresponding to b = 1.0). The differences between nonlinear and linear variants are not great for age in the range of 0–30 years and practically coincide for the age of 30–32 years. It is not surprising because the overwhelming majority of pipes are of the age of 31–32 years (Table 1). However, the great differences for the age of 38–41 years are quite surprising, because the total length of pipes with this age is considerable (Table 1). The analysis of the location of pipes with different age, points of measurement and velocities showed that the pressure was measured in one point only in the oldest part of the WDS and the velocities in the area were very low (0.05 m/ s) even at the highest consumption. That is why both curves (1 and 2 in Fig. 5) give practically the same values of the OF and the surface of the OF has such a shape with thalweg (see Fig. 4). It means that the set of measurements is not sufficient for the evaluation of the exponent with a fair degree of confidence. The points of measurement must be selected by more complicated techniques than in the case of uniform distribution. Therefore, a decision was made to use the linear case of Eq. (2) (exponent b was assumed to be equal to 1.0), to find emax for this variant and use the variant as
1)
14
2)
12 10 8 6 4 2 0
0
10
20 Age, years
30
40
Fig. 5. Dependence of pipe roughness on the age of the pipe with different sets of parameters: (1) emax = 15.3, b = 0.5, OF = 6.80 and (2) emax = 10.3, b = 1.0, OF = 6.81.
the result of calibration. The resulting emax obtained by the LMA with the fixed exponent value b = 1.0 is equal to 10.4 with a standard error of 0.49. Thus, the analysis showed that the visualization of response surface is very useful to understand the results obtained by the LMA. The preparation of such a surface does not take much computer time on the coarse grid (e.g. a grid containing 21 values for each parameter needs 441 calculations by the Epanet2 and takes 4 min on the computer with CPU 3.4 GHz). As was mentioned above, the measurements were divided into two parts. The first part was used to calibrate the model. The second part was used for validation, i.e. the parameters obtained by calibration in the first part of series were used in the calculations for the second part of series, to verify model performance on the independent data. The calculations showed that the sum of squares of differences between the measured and the calculated values was almost the same for the validation of data (the value of the OF was 6.82 for the validation data, as compared to the value 6.80 for the part of data used for calibration). The comparison of the observed and the simulated pressures (validation data) is presented in Fig. 6. The points corresponding to the minimal demand are lower than the other points because at night the pressure decreased in the system.
70
Pressure simulated, m
60
Mean demand Maximal demand
50
Minimal demand
40 30 20 10 0 0
10
20
30
40
50
60
70
Pressure observed, m
Fig. 4. Surface of the objective function on the grid 20 20 points.
Fig. 6. Simulated pressure vs. observed pressure for mean, maximal and minimal demands.
664
T. Koppel, A. Vassiljev / Advances in Engineering Software 40 (2009) 659–664
4. Conclusions The Levenberg–Marquardt algorithm requires a low limit for the step along parameters to calculate partial derivatives in order to work on the operational WDS model containing thousands of pipes. The approximation of the dependence of pipe roughness on age by the function with a low number of parameters is easier than pipe grouping according to age. The analysis showed that the response surface is very useful to understand the results. The combination of the Levenberg–Marquardt algorithm with the response surface obtained by the trial-and-error method facilitates calibration significantly.
Acknowledgements This work was financially supported by Target Financing of Estonia (grant SF0140072s08) at Tallinn University of Technology and by the Estonian Science Foundation (grant G7646). References [1] Vassiljev A, Koppel T. Calibration of the model of an operational water distribution system. In: Topping BHV, Montero G, Montenegro R, editors. Proceedings of the fifth international conference on engineering computational technology. Stirlingshire, United Kingdom: Civil-Comp Press; 2006 [paper 155]. [2] Lansey KE, El-Shorbagy W, Ahmed I, Araujo J, Haan CT. Calibration assessment and data collection for water distribution networks. J Hydraul Eng 2001;127:270–9. [3] Greco M, Giudice GD. New approach to water distribution network calibration. J Hydraul Eng 1999;125(8):849–54. [4] Walski TM. Model calibration data: the good, the bad, and the useless. J AWWA 2000;92(1):94–9. [5] Mallick KN, Ahmed I, Tickle KS, Lansey KE. Determining pipe groupings for water distribution networks. J Water Resour Plan Manage, ASCE 2002;128(2):130–9.
[6] Ormsbee LE. Implicit network calibration. J Water Resour Plan Manage, ASCE 1989;115(2):243–57. [7] Lansey KE, Basnet C. Parameter estimation for water distribution networks. J Water Resour Plan Manage, ASCE 1991;117(1):126–44. [8] Liggett JA, Chen LC. Inverse transient analysis in pipe networks. J Hydraul Eng, ASCE 1994;120(8):934. [9] Vitkovsky JP, Simpson AR, Lambert MF. Leak detection and calibration using transients and genetic algorithms. J Water Resour Plan Manage, ASCE 2000;126(4):262. [10] Kapelan ZS, Savic DA, Walters GA. A hybrid inverse transient model for leakage detection and roughness calibration in pipe networks. J Hydraul Res 2003;41(5):481–92. [11] Methods Haestad, Walski TM, Chase DV, Savic DA, Grayman WM, Beckwith S, et al. Advanced water distribution modeling and management. Waterbury, USA: Haestead Press; 2003. [12] Dennis Jr JE, Schnabel Robert B. Numerical methods for unconstrained optimization and nonlinear equations. Philadelphia: SIAM; 1996. [13] Reddy PVN, Sridharan K, Rao PV. WLS method for parameter estimation in water distribution networks. J Water Resour Plan Manage, ASCE 1996;122(3):157. [14] Echávez G. Increase in losses coefficient with age for small diameter pipes. J Hydraul Eng 1997;2:157–9. [15] Rossman LA. EPANET 2 User’s Manual. Water Supply and Water Resources Division, National Risk Management Research Laboratory, Cincinnati, OH 45268, September 2002. [16] Garcia VJ, Cabrera E, Cabrera E, Jr. The minimum night flow method revisited. In: 8th annual water distribution systems analysis symposium, Cincinnati, Ohio, USA, 27–30 August 2006. [17] Koppel T, Kändler T, Vassiljev A. Calibration of the water distribution network model. In: Brebbia CA, Anagnostopoulos P, Katsifarakis K, Cheng H-D, editors. Proceedings of the water resources management conference, Halkidiki, Greece, September 24–26, 2001. WIT Press: Southampton, Boston; 2001. p. 207–17. [18] Meier RW, Barkdoll BD. Sampling. Design for network model calibration using genetic algorithms. J Water Resour Plan Manage, ASCE 2000;126(4):245–50. [19] Hamming RW. Numerical methods for scientists. 2nd ed. New York: McGrawHill Book Co.; 1973. [20] Gill E, Murray W, Wright M. Practical optimization. New York: Academic Press; 1981. [21] Koppel T, Vassiljev A. Experience in modelling of water supply systems. Congress Publications of ECWATECH-2004, 1–4 June, 2004. [22] Deepan B. Application of optimization methods to calibration of water distribution systems. In: World water and environmental resources congress, June 27–July 1, 2004.