Hydraulic Reliability Assessment and Optimal Rehabilitation/Upgrading Schedule for Water Distribution Systems ()
1. Introduction
Water distribution systems (WDSs) require huge investments in their construction and maintenance. For this reason, it is an important need to improve their efficiency by ways of minimizing their costs and maximizing the benefits. In the last several decades, a significant number of optimization methods have been developed using linear programming, dynamic programming, enumeration techniques, heuristic methods, and evolutionary techniques [1] . However, water distribution systems are so complex that most of the researches only address partial issues. For example, the previous optimal maintenance schedule studies focused on either hydraulic capacity or mechanical performance failure by minimizing costs rather than the combination of them. The proposed study will develop such an innovative tool that considers the deterioration over time of both structural integrity and hydraulic capacity of every pipe, and optimizes the Rehabilitation & Upgrading Schedule (RUS) for an entire WDS.
Hydraulic reliability assessment is used to study the deterioration over time of structural integrity and hydraulic capacity of WDSs; in other words, it is to study the mechanical failure and hydraulic failure over time. However, there are no universally accepted definitions for hydraulic reliability [2] . Here, we define hydraulic reliability as “the ability of a distribution system to meet the demands that are placed on it, where demands are specified in terms of: (i) the flow to be supplied, and (ii) the range of pressure at which these flow rates must be provided [3] .” Hydraulic failures occur when the demand nodes do not receive sufficient flow and/or adequate pressure, due to old pipes with low roughness arising from corrosion and deposition, increases of demand and pressure head requirements, inadequate pipe sizes, insufficient pumping or/and storage capabilities, etc. In addition, due to the fact that the demand is spatially and temporally distributed, hydraulic failures at critical locations in the distribution system may be more important than the average system reliability [4] . Mechanical failures involve pipe breakage, pump failure, power outages, control valve failure, etc. Mechanical failure of the infrastructure components is also integrally linked with hydraulic failure. For instance, pipe breakage is a mechanical failure, and the occurrence of a pipe break will eventually result in insufficient flow rate and/or inadequate pressure head on the demand nodes.
In the previous researches, two methodologies were developed to compute the hydraulic reliability. One is the Monte-Carlo simulation (MCS) technique, which can bring more accurate reliability estimation [5] [6] , but requires a large number of trials to ensure accuracy thereby resulting in a computationally inefficient problem. Another is a linear probabilistic hydraulic model, called the first-order reliability method [7] [8] , which can overcome the computationally inefficient problem in MCS, but has a lower accuracy. Due to the rapid development of computer technology, the MCS has demonstrated its superiority in the hydraulic reliability computation. In this paper, the MCS is chosen to calculate the probabilistic problem of hydraulic reliability. A free simulation software EPANET [9] is incorporated with this MCS to evaluate the nodal & system reliability of a water distribution system. It’s worth noting that the hydraulic reliability assessment here only refers to the hydraulic failure measurement but the mechanical failure measurement.
The RUS optimization model is based on its hydraulic reliability assessment. It is subdivided into two stages. The first stage is the planning of upgrading and paralleling pipes due to hydraulic failure potential (such as demand and pressure); the second stage is the schedule of repairing and replacing pipes due to mechanical failure potential (such as pipe breakage, pump, or aging of pipe). The first stage is domination due to hydraulic failures that directly result in insufficient flow and/or inadequate pressure. On the other hand, mechanical failures only indirectly influence water supply on the demand nodes. For example, the occurrence of one pipe break of a paralleling pipe system may only slightly lower the water flow rate and/or water pressure, but water supply may still be sufficient.
Previously, many research papers used various methodologies to study the optimal maintenance schedule to the pipe networks focusing on either the hydraulic capacity or the mechanical performance failure. Some of the earlier studies used linear programming [10] while the later studies utilized nonlinear programming [8] [11] to study the hydraulic failure for the optimal planning. Many of the recent literature applied genetic algorithms for the determination of the lowest repairing or replacing costs [12] [13] . Some researches had shown several advantages over more traditional optimization methods [14] [15] [16] . In the mechanical failure area, Shamir and Howard applied regression analysis to obtain a relationship for the breakage rate of a pipe as a function of time [17] . Deb et al. suggested a failure probabilistic model to estimate miles of pipes to be replaced on an annual basis [18] . Ascher and Feingold gave a threshold break rate function by using a statistical reliability model [19] .
This paper is to develop such innovative approach, with considers both hydraulic failure and mechanical performance failure, to dynamically assess hydraulic reliability and then continue to find a best option for the long-term RUS of a water distribution system. This approach is conducted on a simulation model of water distribution network in a computer by two universal codes, namely the hydraulic reliability code and the optimal RUS code, and it eventually will be applied in the real water distribution networks planning. The applicability of this approach is verified in an example of a benchmark water distribution network.
2. Methodology
2.1. Hydraulic Reliability Assessment
Hydraulic failures occur when the demand nodes do not receive sufficient flow and/or adequate pressure, due to old pipes with low roughness arising from corrosion and deposition, increases of demand and pressure head requirements, inadequate pipe sizes, insufficient pumping or/and storage capabilities, etc. Mechanical failures involve pipe breakage, pump failure, power outages, control valve failure, etc. Mechanical failure of the infrastructure components is also integrally linked with hydraulic failure. For instance, pipe breakage is a mechanical failure, and the occurrence of a pipe break will eventually result in insufficient flow rate and/or inadequate pressure head on the demand nodes. Mathematically, the system hydraulic reliability equals 1 minus the system failure which shown in Figure 1. It notes that there has not a measurement method to directly identify whether a failure is hydraulic or mechanical.
2.1.1. Nodal Hydraulic Failure Probability
Hydraulic failure at a given node occurs when the supplied flow or pressure is less than
the required minimum flow or the required minimum pressure at that node. Considering we will use EPANET for simulation that the hydraulic simulator always satisfies water demand (because of applying the nodal mass balance equation) but not pressure head, this approach automatically assumes that water demand is satisfied. The nodal mass balance equation can be expressed as:
(1)
where Di is demand at node i; Hi is head at node i; Hj is head at node j; N is total number of the demand nodes with unknown heads; 
is set of nodes connected directly to node i; 
is a nonlinear function relating the hydraulic loss and flow rate in the pipe connecting node i to node j. The nonlinear function used here is the Hazen- Williams or the Colebrook-White equation.
In the above equation, nodal pressure head is an unknown decision variable, and water demands are state variables. This equation can be expressed in a more compact form:
(2)
where Fi(*) is a vector of functions representing the mass balance at each node; Hi is a vector of supplied pressure heads at each node; 
is a vector of state variables and k is the total number of variables. In the study of the hydraulic performance failure, we define water demands Di, pipe coefficients Ci, and tank/storage levels Li as known state variables. Thus, the vector of state variables in this study can be expressed as
; T represents the matrix transpose operation. Moreover, due to uncertainty, those known state variables will be treated as random values. Each random variable is expressed as a probability distribution, and then a random number generator is used to generate the values of Di for each node, the values of Ci for each pipe, and the values of Li for each tank/storage. Incorporation with EPANET, Monte Carlo simulation (MCS) technique can be applied to produce the random values of pressure head.
EPANET: EPANET is free public domain software developed by USEPA. It performs extended period simulation of hydraulic behavior within pressurized pipe networks [9] . A network consists of pipes, nodes, pumps, valves and storage tanks or reservoirs. EPANET tracks the flow water in each pipe, the pressure at each node, the level of water in each tank during a simulation period. Running under Windows, EPANET provides an integrated environment for editing network input data, running hydraulic and water quality simulations, and viewing the results in a variety of formats. In this study, water demands, pipe coefficients, tank/storage levels are inputting data during the simulation; pressure heads are outputs. We assume that pumps and valves are in good operation condition. The failure of pumps and valves is considered in the area of mechanical failure.
Monte Carlo Simulation: Monte Carlo simulation involves repeated generation of pseudovalues for the modeling inputs, drawn from known probability distributions within the ranges of possible values. In this study, this process will be completed by computer, where three random variables (water demand Di, pipe coefficients Ci, and tank/storage level Li) are produced synchronously at each time. The generated pseudovalues are then used as inputs for EPANET model to produce the supplied pressure head values of Hi. The Monte Carlo simulation implementation in a node i includes steps: (1) development of the representative probability distribution functions for selected model input parameters; (2) generation of pseudovalues for each of the selected input parameters from the distribution developed in the previous step; (3) implementation of EPANET model with the pseudovalues to generate a pressure head value of Hi; (4) comparison of the calculated value of Hi with the minimum allowable head
; (5) repeated application of the previous two steps; (6) presentation of the results, and analysis of the results for decision making.
In the above step (4) of the MCS procedure, we defined 
as the minimum allowable head at a demand node. Thus, the probability of hydraulic performance failure at a demand node can be expressed as
(3)
where 
is the probability density function of the supplied pressure head; DiS is the supplied flow rate at the node i, and DiD is the water demand at the node i.
2.1.2. Mechanical Failure Probability
Mechanical failure involves pipe breakages and electric component failures (such as pump failure, power outages, control valve failure, etc). There are many factors that may cause pipe breaks, for instance, aging, physical bending stress, chemical and biological corrosion, etc. Among these factors, aging, material type, dimension, and bedding quality are important factors in predicting pipe breakage. Here we demonstrate two exponential break rate equations. Shamir and Howard’s exponential break rate model describes an exponential relationship between aging and breakage [17] .
(4)
where Aj is the growth rate coefficient of the pipe j; 
is the number of breaks per 1000 ft. length of pipe j in year t; t is time in years; t0 is base year for the analysis (pipe installation year, or the first year for which data are available); M is total number of pipes.
In real situation, a primary cause of unreliability is the removal of water lines or pumps from service primarily due to breaks or power outages. A probability of an electromechanical component 
(k = 1, 2, …, U; where U is total number of electric components), such as a pump or a valve, being out of service can be included as another random variable that supported by historical data.
2.1.3. System Hydraulic Reliability Assessment (See Figure 2)
Different water demand has different importance. For example, a hospital water demand has the highest importance, and the water demand of a park area may have a lower importance. We weight failure probabilities to reflect their different importance
![]()
Figure 2. Procedure of system hydraulic reliability assessment.
impacts. Then, the system failure probability 
can be defined as the maximum individual failure probability in the system
(5)
where 
are weighting coefficients; ![]()
is the probability of hydraulic performance failure at the node i; ![]()
is the failure probability of pipe j; ![]()
is the failure probability of electromechanical component k; ![]()
is the fire flow failure probability.
The system hydraulic reliability ![]()
is given by
![]()
(6)
2.2. Optimal Rehabilitation & Upgrading Schedule
The long term planning of WDSs is an important issue for municipal agents. These agents pay special attention to improve the efficiency of their planning by the ways of minimizing the costs and maximizing the benefits. However, water distribution systems are so complex that most of the researches only address partial issues. For example, the previous optimal maintenance schedule studies only focused on either hydraulic failure or mechanical failure by minimizing costs but the combination of them. This study will develop such an innovative model that considering the deterioration over time of both structural integrity and hydraulic capacity of every pipe, and optimizing the RUS for a whole WDS.
The model identifies the network maintenance into the rehabilitation and the upgrading. It is subdivided into two stages (see Figure 3). The first stage is the planning of
![]()
Figure 3. Flowchart of two-stage method.
upgrading and paralleling pipes due to hydraulic failure potential (such as demand and pressure); the second stage is the schedule of repairing and replacing pipes due to mechanical failure potential (such as pipe breakage, pump, or aging of pipe). The first stage is domination due to hydraulic failures that directly result in insufficient flow and/or inadequate pressure. On the other hand, mechanical failures only indirectly influence water supply on the demand nodes. For example, the occurrence of one pipe break of a paralleling pipe system may only slightly lower the water flow rate and/or water pressure, but water supply may still be sufficient.
It’s worth noting that we use “potential” to describe hydraulic failures and mechanical failures because we refer to those failures in future. Both prediction models of hydraulic failure and mechanical failure are based on the hydraulic reliability assessment. The hydraulic failure probability is assessed by the Monte Carlo Simulation, and the pipe breakage probability is computed by the Shamir and Howard’s exponential break rate model (shown in Equation (4)). Because mechanical failure probability is obtained in terms of historical data, the optimal upgrading results will not affect the pipe breakage probabilities except for the pipes that had been upgraded in the first stage.
2.2.1. Rehabilitation Optimization Model
Although the pipe repair cost is much lower than the replacement fee, frequent breaks will result in the cumulative repair expense exceeding the replacement cost. Therefore, we need a methodology to find the critical (or threshold) break rate to determine the timing of repair or replacement [19] .
If we assume that the pipe will be replaced at the time of the nth break (it also implies that for the previous (n − 1) breaks only repairs have been performed), we can write the present worth of the total cost of the pipe as
![]()
(7)
where ![]()
is the discount rate; ![]()
is the time of the ith break measured from the installation year; ![]()
is the repair cost of the ith break; ![]()
is the replacement cost at time![]()
; ![]()
is the present worth.
For the total cost ![]()
at time ![]()
to be a minimum, it must satisfy the condition
![]()
(8)
However, the critical check is to determine the first instance when the condition ![]()
holds true. Solving for![]()
, we have
![]()
(9)
Recognizing ![]()
is the time between the nth and (n + 1)th breaks or time interval for the occurrence of one break, we obtain the threshold break rate, Brkth, as the inverse of ![]()
(where![]()
). The threshold break rate is defined as the break rate between subsequent breaks.
![]()
(10)
Now, from the observed data for any given pipe we can derive a current break rate. Whenever the current break rate is equal to or more than Brkth, the pipe should be replaced. It is worth noting that the optimal rehabilitation schedule information will be transferred to the upgrading optimization process (shown in Figure 4).
2.2.2. Upgrading Optimization Model
The hydraulic reliability can indicate the hydraulic performance information for each demand node. If the hydraulic failure probability on a node is lower than the allowable probability, the water provider should replace the inlet pipe with a larger diameter pipe or build a paralleling pipe. Then, the next few questions for the decision maker are: 1). Replace the existing pipeline by a bigger diameter pipeline or build a paralleling pipeline? 2) What size of pipe should be applied in this replacement or paralleling? 3) The duration of the new inlet pipe that can transfer sufficient water?
The design of water distribution systems is often viewed as a least-cost optimization problem in the previous literatures. The pipe diameters were considered as decision variables. Obviously other parameters should be considered as objectives in the optimization process, such as reliability, redundancy and water quality [20] . However, problems with quantifying these objectives for use within optimization design models kept researchers concentrating on the single, least-cost objective. Therefore, the overall objective in the upgrading optimization in this paper is to minimize the total cost. The following formulation describes this objective:
![]()
Figure 4. Example of distribution system.
![]()
(11)
where ![]()
is the cost of pipe j with diameter Rj and length Lj; m is the number of pipes that must be upgraded. The objective is subject to the following constraints:
![]()
(12)
![]()
(13)
![]()
(14)
![]()
(15)
where Equation (12) refers to conservation of flow constraints and energy equations (loop equations); Equation (13) is head boundary constraints; Equation (14) is designed water demand boundary constraints; Equation (15) represents general constraints including financial constraints.
The decision on either replacement or paralleling pipeline is depended on cost, geographical situation, and the existing pipe condition. Some pipeline geographical conditions are not suitable for the laying of parallel pipelines that the only choice is to replace the existing pipelines by larger diameter pipes. However, the geographical condition consideration of the pipe replacement planning is out of the scope of this study. In addition, we have to consider the duration of the upgraded pipeline. If one pipe needed upgrade because low hydraulic capability, simultaneously it was aging pipe that would be replaced in the next few years, it was better to directly be replaced by a new large diameter pipe rather than to build a parallel pipeline.
We assume the length of the term is l, the objective formulation is modified as follows
![]()
(16)
where l is the length of the planning period; and k (k = 1, 2, … l) is the year index of pipeline upgrading. The total cost includes the following components:
![]()
(17)
![]()
represents the pipeline costs including installation, paralleling, replacement and repair costs; ![]()
is the cost of setting up the construction plant and machinery.
2.2.3. Optimization Procedure
Pipe Roughness and Water Demand: The pipe roughness and the water demand are two model input parameters that vary over time. The pipe roughness increases over time at a rate that varies according to the pipe type, water quality, and operation and maintenance practices. To model the effect of aging on the carrying capacity of pipes, the equation of Sharp and Walsli [21] was used
![]()
(18)
where C(t) is the Hazen-Williams coefficient in year t; eo is the roughness at the time of installation (mm); a is the roughness growth rate (mm/year); and R is the pipe diameter (mm).
The function developed by Dandy et al. [22] is used to predict nodal water demand.
![]()
(19)
where D0j is the base year water demand for the node j and DGR is the annual percentage rate of increase in the base demand. DGR was assumed to follow demographic data patterns closely, and so it was taken as the population growth rate with typical values of about 2 to 5; t is the time in years; Pt and P0 are the price per unit volume of water for year t and the base year, respectively; PREL is the price elasticity of demand. A typical range of PREL values is −0.2 to −0.5. Generally it is advisable to increase the price of water in the last 1 or 2 years of the design period for the best results in order to delay the need to upgrade the network. When there is adequate capacity after upgrading, prices can be reduced in order to utilize the system capacity as fully as possible and to maximize economic efficiency. Other demand management options could be used, and self-evidently, the joint effect of all the demand management techniques deployed should be considered.
Genetic Algorithms: We have chosen the genetic algorithms (GA) as the optimization method. GA is a search technique used in computing to find exact or approximate solutions to optimization and search problems. It is based on the mechanics of natural selection and natural genetics that using biological techniques including inheritance, mutation, selection, and crossover [23] . GA provides a robust and efficient way to search complex parameter spaces for ever better solutions to an optimization problem [23] . Although there are no guarantees that a GA will actually attain the optimum, it generally finds one or more extremely good solutions with relatively little computational effort.
The GA code can also be written by MATLAB. Moreover, EPANET simulation will be run during the GA optimization process. The following steps give the details of the optimization procedure for a repair schedule.
Step 1. Use Equation (4) to calculate each pipeline’s break rate B(i, j); and use Equation (10) to compute each pipeline’s Brkth and then obtain the corresponding replacement year Yr(j), where i represents the time of year, and j is the pipe index.
Step 2. Set i = 0;
Step 3. i = i + 1;
Step 4. Run the hydraulic reliability model and calculate nodal hydraulic reliability R(i, j) for all studied nodes in a WDS in the ith year;
Step 5. Compare the predicted hydraulic reliability with the minimum allowable hydraulic reliability Rm(i, j) at the ith year, to find the pipes to be upgraded P(i, j);
Step 6. If the replacement year Yr(j) of the upgraded pipe is equal to or more than the planning term N, then, no paralleling is considered;
Step 7. Run Optimization Code (GA) to find optimal diameter pipe P(i, j, m) (where m is the index of pipe diameter, both replacement and paralleling use this diameter index), with the objective of minimizing total cost, and subject to the constraint that the R(i, j) should be greater than the maximum of Rm(i, j) during the whole planning period N. The computation of R(i, j), of course, needs to run the hydraulic reliability prediction model;
Step 8. If i < N, then go to step 3;
Step 9. Calculate the total cost and record the optimal rehabilitation and upgrading schedule.
3. Example Application
Figure 4 shows a benchmark hypothetical water distribution system [24] that consists of 16 pipes, 10 demand nodes, two tanks and one reservoir. Tables 1-4 demonstrate the information of system components.
Due to we use the MCS to compute the nodal and system reliability, the iterative process of random number generation of D, C, and L, and hydraulic simulation must be repeated a large number of times. According to the previous studies [6] [8] , the optimal number of iterations is between 500 times to 3000 times. Because the results had tiny different after 2000 iterations, we use 2000 iterations in this study.
There are several probability distributions can achieve random number generation of D, C, and L, such as normal, log-normal, uniform, Pearson and Weibull distributions. According to the literature, normal distribution is the most popular probability distribution that used in hydraulic studies [13] [14] [25] [26] . We also use normal distribution for random number generation of D, C, and L.
The Shamir and Howard’s exponential break rate model (shown in Equation (4)) is used to compute the pipe breakage probability. The expected number of failures per year per unit length of pipe Aj is uniformly assumed to be 0.122. The mean values of
normal distribution of D, C, and L are obtained from Table 1 to Table 4. The coefficients of variation for D and L are assumed to be 0.2. According to the conclusion of Bao and Mays [6] that system reliability is somewhat insensitive to the type of probability distribution of pipe roughness when the coefficient of variation C is less than 0.4, we assume the coefficient of variation C is 0.4.
The results of the nodal hydraulic performance failure probability, the pipe breakage probability, and the system hydraulic reliability are shown in Table 5. From the results, we found i). The system reliability shows a great improvement when a reasonable weighing factor is used. ii). The system reliability is 0.892, representing a sample distribution of an unsafe system (safe system reliability should be greater than 0.999). iii). the nodal hydraulic performance failure probability is almost one order of magnitude lower than the pipe breakage probability.
![]()
Table 5. Results of reliability assessment.
We still use the water distribution system illustrated in Figure 4 to demonstrate the optimization method for repair schedule. Assuming that this is a 15-year water distribution system maintenance plan that with the following cost data (price on the first year). It includes the total cost of pipe, installation, construction and labor fees (shown in Table 6).
If the discount rate is 0.06 (shown in Table 7), according to Equation (10) and the cost data in Table 6, we obtain the threshold break rate for each pipe. Moreover, using Equation (4), we can calculate the replacement timing for each pipe.
Finally, after running the optimization model in terms of those data, an optimal repair schedule is shown in Table 8. It indicates that the number of pipe break grows as time increases until the pipe replacement conducts. The first year has 5 pipe breaks, the number of pipe break increases to 8 in the fifth year, and only two pipe breaks happen in the last year of the schedule term due to the pipe upgrading. The first pipe replacement
![]()
Table 7. Threshold break rate and replacement timing (discount rate = 0.06).
![]()
Table 8. The optimal repair schedule.
occurs in the ninth year that pipe P9 is replaced by an 8-inch diameter pipe. Only one paralleling pipe is built in the fifteenth year of the schedule term. The total cost of this schedule is $471,933.
4. Conclusion
The paper successfully developed an innovative approach to optimize the rehabilitation & upgrading schedule for water distribution systems based on the stochastic hydraulic reliability assessment. The approach considered the uncertainty of hydraulic parameters and treated them as random numbers following the normal distribution. The hydraulic probability was computed by using the Monte Carlo simulation technique and incorporated EPANET simulation software in MATLAB. This approach subdivided the high complexity WDS optimization problem into two stages―one stage for hydraulic performance optimization and another for mechanical performance optimization. The optimizations of two stages were performed separately but exchanged information with each other. Simultaneously, the threshold break rate model and the genetic algorithm participated in this optimization process. Two practical codes, the hydraulic probability assessment and the RUS schedule, were developed based on this approach and applied for a benchmark water distribution system. The results demonstrated the applicability of this approach.