Iterative hydraulic interval state estimation for water distribution 1 networks

17 State estimation of hydraulics (i.e. pressure and flows) in water distribution networks is 18 an important tool for efficient and resilient operation. However, hydraulic state estimation 19 is a challenging task in practice due to the scarcity of measurements and the presence of sev20 eral modeling uncertainties. Standard state estimation techniques may produce unreliable 21 estimates with no information of the estimation error magnitude, especially when historical 22 data are used in place of missing measurements. In this paper, we propose a comprehensive 23 1 Vrachimis et al., August 9, 2017 methodology for generating hydraulic state bounding estimates by considering both mea24 surement and parametric uncertainties. The methodology is based on solving the nonlinear 25 interval hydraulic equations using bounding linearization, a technique that restricts the non26 linearities within a convex set, thus converting the problem in a form which is solvable using 27 linear optimization. An iterative procedure improves the bounding linearization, converging 28 to the tightest possible bounds. Simulation results demonstrate that the proposed method29 ology produces tight state bounds that can replace Monte-Carlo simulations. 30


INTRODUCTION
The water industry is being modernized with the installation of sensors for monitoring Water Distribution Networks (WDN) and computer systems to process these data.Integrated platforms are already being developed that algorithmically combine real-time sensor measurements from Supervisory Control and Data Acquisition (SCADA) systems, geographic information systems (GIS), and hydraulic models to provide useful information to the operators.A state estimation algorithm infers the complete system state, such as water flows in pipes, consumer water demands, pressures at nodes and tank levels, using the available measurement set and network equations.A complete view of the network state supports the decision-making process and enables the efficient operation of these systems, improves customer service and enables the early detection of emergency events, thus minimizing their impact.Examples of the use of state estimation in real systems include the use for online burst detection (Okeya et al. 2014) as well as for online modelling (Machell et al. 2010) and control of WDN (Rao and Salomons 2007).
Standard state estimation techniques require a measurement set that makes the system observable, i.e. the sensor number and locale ensure that the system state can be calculated (Bargiela 1985;Nagar and Powell 2000;Díaz et al. 2016).Additionally, the statistical characterization of sensor measurement error is needed to give more weight to measurements originating from more accurate sensors.Then, using a mathematical model of the network, a state estimation algorithm can infer the system state.Many approaches have been proposed to solve the state estimation problem for water systems, such as Kalman Filtering and Weighted Least Squares (WLS), with the latter being the most widely used and varied.The above methods produce a point in state-space and are referred to as point state estimation (Powell et al. 1988;Andersen et al. 2001;Kang and Lansey 2009).
State estimation in WDN is a challenging task, mainly due to the scarcity of measurements.In contrast to other large scale engineering systems (e.g.power systems) where the number of measurements guarantees observability, in WDN observability is almost never achievable with the available measurements.Some parts of the WDN may be widely monitored, such as the transport network, however even a single sensor failure could make these parts unobservable (Vrachimis et al. 2016).The large area covered by WDN and the large number of system states is one of the main reasons that an impractical number of sensors must be installed to guarantee observability.A common practice to reduce the complexity, is to skeletonize the network by treating a group of consumers as a single demand point.It is then possible to use pseudo-measurements, which are demand estimates determined from population densities and historical data, to complement the missing measurements (Hutton et al. 2014).Recent advances in water demands research have also made possible the higher resolution modeling of water demands, thus reducing the need for skeletonization of networks and increasing the accuracy of state estimation (Avni et al. 2015).
The use of pseudo-measurements may introduce new problems to the state estimation process, as they are highly uncertain and the resulting estimates may deviate significantly from the real system state.This in turn could affect other algorithms which rely on stateestimation, such as feedback control or fault-diagnosis.Efforts have been made to characterize the uncertainty of pseudo-measurements (Bargiela and Hainsworth 1989), but it is improbable that a statistical characterization will be available.Thus, standard state estimation techniques such as WLS may not be capable of producing a reliable measure of the estimation error.Consequently, researchers have tried to combine online estimation of demands with state estimation, in order for the latter to be more accurate (Preis et al. 2011).
Another significant source of uncertainty which complicates WDN state estimation is modelling and parameter uncertainty.Recent works provide explicit expressions for the sensitivity of the state estimation problem to these uncertainties (Díaz et al. 2018).Typically, the network topology is assumed known, especially after the process of skeletonization which simplifies the network graph.However, even when the topology is known, pipe parameters such as length and diameter are rarely known accurately and estimates are used in place.This is especially true for pipe roughness coefficients, which along with length and diameter, are used to calculate the headloss across pipes.This is why, even with an observable sensor configuration, model calibration is required a priori or during state estimation for the latter to produce feasible solutions (Gao 2017).Model calibration can be considered as the inverse problem of state estimation, during which estimates of the unknown model parameters are calculated based on measurements (Kapelan et al. 2003;Savic et al. 2009).But even after calibration, it is possible that the calculated parameter set satisfies the constraints imposed by measurements, but deviates from the true parameter set.
Considering the many unknowns and uncertainties in WDN state-estimation, it is evident that accurate state-estimates are difficult to be generated without some kind of trade-off.A practical approach for state estimation in the presence of demand and modeling uncertainty, is interval state estimation (Bargiela et al. 2003;Langowski and Brdys 2007).This approach models the uncertainties on input data as intervals, defined by lower and upper bounds.
Then, considering this bounded uncertainty, interval state estimation provides lower and upper bounds on the state estimates, in contrast to point state estimation methods which only provide a single point.Providing a range of values for each state, is often more useful to an operator than providing point estimates which give no indication on their proximity to the true state value.Additionally, having reliable interval state estimation is essential in many methodologies related to event and fault detection such as leakage detection (Pérez et al. 2009), water contamination detection (Eliades et al. 2015) and sensor fault detection (Vrachimis et al. 2015).
The use of bounds for the representation of measurement uncertainty and the subsequent calculation of state estimate bounds was introduced in (Bargiela and Hainsworth 1989).This idea was further developed in (Brdys and Chen 1994) as the set-bounded state estimation problem.The process of calculating bounds for state estimates caused by measurement uncertainty is also referred to as Confidence Limit Analysis which can be solved using different approaches, including Neural Networks (Gabrys and Bargiela 1997), the Error Maximization method (Arsene et al. 2011), the Ellipsoid method and Linear Programming (Bargiela et al. 2003).All these approaches assume a known network model which can be linearized in order to solve the non-linear equations that characterize WDN and provide state bounds based on measurement uncertainty.Few methodologies can guarantee the inclusion of the true state in the bounds based on given uncertainty, while the effect of modeling uncertainty is not considered.Another approach that could incorporate modeling uncertainty is the use of Monte-Carlo Simulations (MCS), where state bound estimates are obtained by randomly generating and evaluating a large number of model parameter sets or realizations (Pasha and Lansey 2010).This approach requires a sufficiently large number of simulations, and even then some cases may not be covered, leading to underestimation of the range of the true state bounds.
In this work we propose a new interval hydraulic state-estimation approach for WDN that considers the combined effect of bounded measurement and modeling uncertainties.The proposed methodology calculates the bounds on state estimates using the nonlinear form of the network equations, by also modeling pressure-dependent demands and background leakages.The nonlinear modeling guarantees that when accurate uncertainty bounds are provided, the bounds on state estimates will include the true system state.This is achieved using bounding linearization, a technique which restricts the nonlinearities within a convex set, thus converting the hydraulic equations in a form where the minimum and maximum of each state can be found using linear optimization.Then, an iterative procedure is followed to minimize the distance between upper and lower state bounds, by improving the bounding linearization at each step and converging to the tightest possible bounds.The contributions of this work are: • The consideration of both modeling uncertainty, in the form of uncertain parameters, as well as measurement uncertainty in the interval state estimation problem for WDN.
• The development of a novel algorithm that calculates tight hydraulic bounding estimates based on the considered uncertainties.
• The use of the nonlinear form of the hydraulic equations which also considers pressuredependent demands and background leakages, in order to ensure that the bounding estimates guarantee the inclusion of the true system state if the uncertainties have been accurately represented.
This paper is organized as follows: Section "Problem formulation" formulates the problem of hydraulic state estimation in WDN where the uncertainty on model parameters and measurements is represented by intervals.Section "Iterative Hydraulic Interval State Estimation" presents a methodology to solve this problem based on the Iterative Hydraulic Interval State Estimation (IHISE) algorithm.In Section "Case Studies" this methodology is applied on different benchmark water networks and its performance is assessed.

PROBLEM FORMULATION
The topology of a WDN is modeled by a directed graph denoted as G = (N , L).Let N = {1, • • • , n n } be the set of all nodes, where |N | = n n is the total number of nodes.These represent junctions of pipes, consumer water demand locations, reservoirs and tanks.The unknown quantity associated with nodes is the hydraulic head, indicated by h i .Reservoirs and water tanks that have level sensors installed, can be considered as nodes with known head.We define the set of nodes with unknown head is the number of nodes with unknown head.The set of nodes with known head is defined as is the number of nodes with known head and N = N u ∪ N h .Each node j with unknown head is associated with a water consumer demand at the node location, denoted by q ext,j .
Let L = {1, • • • , n l } be the set of links, where |L| = n l is the total number of links.
These represent network pipes, water pumps and pipe valves, with the last two being the main hydraulic control elements in a water network.We define the set of links that represent pipes as L p = {1, • • • , n p }, where |L p | = n p is the total number of pipes.We also define the set of links that represent pumps as L pu = {n p + 1, • • • , n l }, where |L pu | = n pu is the total number of pumps.The unknown quantity associated with a link i is the water flow, indicated by q i .

Formulation of hydraulic equations
It is a common practice in WDN to receive sensor measurements of flows, pressures or tank water levels at constant time intervals, which typically range from five minutes to one hour.These sensors may also give an average measurement for the elapsed time interval, thus fast changing dynamics (e.g.pressure transients) are neglected.As a result, standard state estimation in WDN is carried out at steady state, with the system state being recalculated when measurements arrive.
In this work we assume that only lower and upper bounds on measurements are available.
The bounds can be derived from real sensor measurements, or from population densities and historical data (pseudo-measurements).The measurement bounds are available at every discrete time step k for all nodal demand outflows and for all tank and reservoir levels.This sensor configuration guarantees the topological observability of the network.Other sensor configurations are also possible, given that they satisfy the topological observability condition, which can be checked using the algorithm in (Díaz et al. 2017).The vector of measured external water demand outflow is indicated by qext (k) ∈ R nu .The known head vector, which results from tank and reservoir level measurements, is indicated by hext (k) ∈ R n l .The unknown state of the network are the water flows in pipes q(k) ∈ R n l and the unknown hydraulic heads at nodes h(k) ∈ R nu .
The state is calculated using a hydraulic model of a WDN, which is a set of equations derived from the laws of: 1) conservation of energy and 2) conservation of mass in the network.In this work we use the pipe formulation of these equations as used by (Todini and Pilati 1987), which has been shown to be robust in computer simulations (Rossman 2000).
The only dynamic component of these equations are the changing tank levels (Boulos et al. 2006).Because tank levels are assumed to be measured, the resulting hydraulic equations are not dynamic, thus the discrete time notation k is omitted.The formulation of these equations follows.

Conservation of energy equations
Energy in WDN is associated with the head at nodes and when water flows through a network link i which connects two nodes, a flow dependent head function f i (q i ) describes the change in head.In the case of pipes, energy is dissipated due to friction of water flowing through the pipe, resulting in head-loss between two connected nodes.Head-loss depends on the water flow through the pipe but also on pipe parameters.Each pipe i ∈ L p is characterized by pipe length l i , pipe diameter d i and pipe roughness coefficient c i .All these quantities are used in the empirical Hazen-Williams (H-W) formula (Boulos et al. 2006) to calculate head-loss.The effect of pipe parameters in this formula is aggregated in the H-W resistance coefficient r i of each pipe, which is a function f r HW : R + × R + × R + → R + of pipe parameters, defined as: The head-loss across pipe i ∈ L p is then calculated using the H-W formula as follows: where ν is a constant exponent associated with the H-W formula and q i is the water flow in pipe i.
Another example of network element are pumps i ∈ L pu which are characterized by a head-flow curve, which is used to relate the flow through the pump to the head-gain across the pump, according to each pump specifications.This is given by: where w 1 , w 2 , w 3 ∈ R are coefficients of the pump head-flow curve, while the minus sign indicates that in the case of pumps there is head-gain instead of head-loss.
The energy equations for all the network links, considering elements modeled by ( 1) and (2), can be written as follows: where: • h j is the unknown head of node j ∈ N u .
• B ∈ R n l ×nu is the incidence flow matrix, indicating the connectivity of nodes with links.Element B ij = +1 if the direction of link i enters node j; element B ij = −1 if the direction of link i leaves from node j; otherwise B ij = 0. Nodes with known head are excluded from this matrix.
• h ext,i is the sum of known (measured) heads that appear in each equation i ∈ L. In vector notation, the known head vector is given by hext ∈ R n l .

Conservation of mass equations
The conservation of mass law for a node j ∈ N u is similar to Kirchhoff's current law in electric circuit analysis and can be summarized as follows: the sum of branch water flows from pipes incident to a node j must be equal to the node's external water demand q ext,j .
A demand-driven modeling approach assumes that the demand at each node is independent of the pressure at that node.However, this analysis is not valid when power outages, fire fighting, pipe breaks or temporarily closed portions of a WDN lead to pressure-deficit conditions.In those cases, the consumers do not receive the requested demand, thus the modeling of demand is no longer valid and a pressure-dependent demand modeling is recommended.
The pressure-demand relationship can be modeled by multiplying the user requested demand q ext,j at node j by the pressure depended function f ext,j (h j ), which is given by (Wagner et al. 1988;Giustolisi and Laucelli 2011;Klise et al. 2017): In ( 4), H req,j is the head above which the consumer can receive the requested demand q ext,j (depends on node elevation), H min,j is the minimum desired head at node j (depends on node elevation) below which the consumer does not receive any water.
Background leakage flows are also present in real WDN and are modeled as an added demand component at nodes.Leakage flows are pressure-depended and are modeled similarly to pressure-driven demands as follows (Giustolisi et al. 2008): where Z j is the elevation of node j, and β j and γ j are leakage parameters depending on pipe deterioration and material.
The conservation of mass equations, considering all the nodes of the network, can be written using the incidence flow matrix as follows: i∈L In vector notation, the requested external water demands for all nodes are given by qext ∈ R nu and the leakage flow at nodes by qleak ∈ R nu .Equations ( 3) and ( 6) define the network state, which are the water flows in pipes and hydraulic heads at nodes, indicated by x = [q h ] .

Measurement and parameter uncertainty
As mentioned in the "Introduction" section, we consider sensor measurements, or pseudomeasurements that are uncertain, with a bounded measurement error.We assume that these are available for all nodal demand outflows, as well as for all tank and reservoir levels, guaranteeing an observable sensor configuration.The network topology is available in the hydraulic model, but pipe parameters are only approximately known.
The uncertainties are modeled in this work using intervals, which is equivalent to a uniform probability distribution.For notational convenience, we adopt the convention of denoting intervals in bold.Let v = [v l , vu ] be a closed interval vector, where vl is the lower bound vector and vu is the upper bound vector, such that: ., n}}, and n is the size of the vector.Note that calculations performed in equations containing intervals require the use of interval arithmetic (Daumas et al. 2009;Kearfott 1996;Moore et al. 2009).
The uncertain requested water demands and the reservoir/tank levels are given by the interval vectors qext = q l ext , q u ext and hext = h l ext , h u ext respectively.Note that, measurement bounds that are derived from nodes with actual sensors, do not require to be accompanied by a pressure-dependent function, so f ext,j (h j ) = 1.
We also consider the uncertainty on the head function f i (q i ).When this function contains uncertain parameters, these will be modelled as intervals defined by a lower and upper bound, and the head function will be indicated in bold as f i (q i ).Uncertainty in pipe parameters is included in the uncertain H-W coefficients r i .These are calculated using uncertain pipe parameters, which are the roughness coefficients c i , diameter d i and length l i .For a certain pipe i, the uncertain H-W coefficient is given by: r i = r l i , r u i .An interval H-W coefficient transforms the pipe headloss function given by (1), into an interval function given by: Similarly, for i corresponding to a pump with an uncertain pump curve, (2) becomes: The uncertainty is considered in the case of leakage modeling, is on the leakage parameters, which can be represented by the intervals β j = β l j , β u j and γ j = γ l j , γ u j .This will result in an interval leakage function q leak,j (h j ).
As a practical note, a calibration pre-step of network topology parameters is recommended, as it will reduce the uncertainty of these parameters in the sense that their boundaries will be more restrictive.As a result, the bounded state estimates calculated by the IHISE algorithm will be less conservative.In this work the parameter boundaries are assumed constant for all time steps, because these parameters vary very slowly over time.The parameters can be updated whenever a calibration procedure takes place for the network.

Problem definition
The problem of solving the hydraulic equations of a WDN when these contain uncertainty in the form of intervals, is reduced to finding the set of all point solutions for the state x = [q h ] , that satisfy the following systems of equations: Problem ( 9) is a Nonlinear Interval System of Equations (NISE) and the set of solutions for x that satisfy (9) may have a rather complex form that needs to be described with nonlinear functions.This is why, in the literature, 'interval solutions' are most often considered, with the aim of finding the Interval Hull (IH) solution, i.e. the smallest interval vector x containing all solutions for the NISE.The literature on finding an OI solution to NISE is limited, but some approaches have been proposed, such as (Kolev 2004b), which uses affine arithmetic to represent the equations and interval linearization to deal with the nonlinearities.However, this approach does not consider interval multiplicative terms in the nonlinear functions, thus cannot be applied to (9).The solution to NISE can also be approached using optimization, as in (Jiang et al. 2008) where the task of solving nonlinear interval number programming problems was investigated.
This method, however, does not ensure that the solution is an OI, thus it is not suitable for use with methodologies such as fault detection which require outer bounds on states.
Good approaches in the literature that provide tight OI solutions exist for Linear Interval Systems of Equations (LISE) and are mainly divided in two categories.The first uses interval arithmetic (Daumas et al. 2009;Moore et al. 2009) to find the solution.Due to the fact that when using interval arithmetic to solve LISE, the solution interval is inherently an overestimation, iterative methods are used to approximate the IH solution, such as the Gauss Elimination method, LU decomposition method and the iterative Jacobi method (Zieniuk et al. 2015).The second category formulates LISE as an Interval Linear Programming problem, where intervals can exist both in the objective function and in the constraints (Chinneck and Ramadan 2000;Huang and Cao 2011).This approach is promising since the formulation of the equations as a Linear Program provides the opportunity to add and manipulate constraints to the problem.Additionally, powerful software that solve Linear Programs efficiently exist, which reduce computation time.

ITERATIVE HYDRAULIC INTERVAL STATE ESTIMATION
In this work we propose an iterative method for finding an OI solution of the NISE in (9), named Iterative Hydraulic Interval State Estimation (IHISE), which closely approximates the IH solution.This method deals with the nonlinearities in (9) using bounding linearization, which encloses the interval nonlinearities in a convex set and converts (9) into a system of linear inequalities.The resulting linear inequalities are then appropriately formulated into a Linear Program and new bounds on the state variables are calculated.An iterative procedure then approximates the IH solution of (9) by using the new bounds on the states to improve the bounding linearization.Initial bounds on state variables can be defined from physical properties of WDN.
The IHISE algorithm is comprised of five main steps: 1. Find initial bounds on the state variables using physical constraints of the system.
2. Use bounding linearization to bound the nonlinearities in a convex set.
3. Formulate Linear Programs (LPs) for each state using the resulting linear inequalities.
4. Solve one maximization and one minimization LP for each state to produce new upper and lower bounds.
5. Iteratively improve the OI solution of ( 9) by repeating steps 2 to 4 until convergence of bounds.
Next, the five steps of IHISE are described in detail.
Step 1: Initial bounds on state vector The first step of the IHISE algorithm is to impose initial bounds on the state vector x = [q h ] which should be an OI solution of ( 9).The initial bounds on the unknown head vector h are given by the interval vector h(0) = hl (0) , hu (0) and are chosen by considering physical properties of the network.The minimum head vector hl (0) is set equal to the elevation of each node and the maximum head vector hu (0) is set equal to the sum of reservoir and pump heads, which is the maximum head that any node in the network can have.
The special structure (9a), in which each equation contains only one flow state q i , allows us to use the initial bounds on heads h (0) to find the initial bounds on the flows.We rewrite (9a) as follows: where y i = y l i , y u i is a known interval derived from the known terms in (10) using interval analysis.The function f i (q i ), when i ∈ L p , is given by ( 7).This function is inclusion isotonic (Moore et al. 2009) meaning that if q 1 i ⊆ q 2 i then f (q 1 i ) ⊆ f (q 2 i ).This property enables the derivation of analytical bounds on the unknown pipe flows, by rearranging (10) with respect to q i , i ∈ L p .
In the case of pumps, f i (q i ) is given by ( 8).This function is not inclusion isotonic in its whole range, but this property holds in the special case when q i > 0 or q i < 0. This implies that the flow direction in pump links needs to be known a priori, which is a valid assumption, since pumps are physically restricted to move water in one direction.Assuming that the flow in pump links is always positive, the bounds on flows q i , i ∈ L pu can be found by rearranging (10) with respect to q i .
The initial bounds on the flow state vector are denoted by q(0) = q l (0) , q u (0) .An analytical derivation of these bounds for pipes and pumps is given in Appendix S1 of the Supplemental Data.
Step 2: Bounding linearization of interval nonlinear terms This step aims at enclosing in a convex set S the nonlinear uncertain functions that exist in Problem (9).This will allow the formulation of a Linear Program.Problem (9) contains three nonlinear uncertain functions: f i (q i ), q ext,j f ext,j (h j ) and q leak,j (h j ), which are all functions of one bounded variable.The bounds on these variables have been calculated Vrachimis et al., August 9, 2017 in Step 1, such that for flow variables q i ∈ q l i , q u i and for head variables h j ∈ h l j , h u j .The goal is to construct convex sets S that include all the points of the uncertain functions in the given range, e.g.f i (q i ) ∈ S, ∀q i ∈ q l i , q u i , i ∈ L.
In this work we design the convex sets S using bounding linearization (Kolev 2004b).
This procedure can be used on any uncertain nonlinear function of one bounded variable and it encloses the function between two lines.For example, given the nonlinear uncertain function f i (q i ) for an interval q i ∈ q l i , q u i , a lower line f l L,i (q i ) = λ l i q i + β l i and an upper line f u L,i (q i ) = λ u i q i + β u i can be designed such that: The goal of the bounding linearization procedure is to define the line parameters to minimize the area of the resulting convex set S. A detailed description on how to obtain linearization bounds for each of the nonlinear functions considered can be found in Appendix S2 of the Supplemental Data.
Step 3: Formulation of Linear Program The bounding linearization of Step 2, yields linear inequality constraints for the interval functions f i (q i ), q ext,j f ext,j (h j ) and q leak,j (h j ).These inequalities can replace these functions in (9) with new variables ζ e,i , ζ p,j and ζ l,j respectively, thus transforming these equations into linear inequalities.Bound inequalities also arise on state variables q i and h j through Step 1.
The Linear Program will then have the following constraints: λ l e,i q i + β l e,i ≤ ζ e,i ≤ λ u e,i q i + β u e,i , i ∈ L (12c) Note that the interval parameters in (9a) are eliminated through the use of their upper and lower bounds to convert them into the inequalities (12a) and (12b).The LP decision variables vector is defined as z = [x , ζ ] ∈ R (2n l +3nu) where x = q , h is the state vector and ζ = [ ζe , ζl , ζl ] is the auxiliary variable vector.Using the constraints (12a) -(12g), two LP problems can be formulated for obtaining lower (LPmin) and upper (LPmax) bounds on each state z i of the vector z: Step 4: Solution of the linear interval system of equations The objective of the optimization problem formulated in the previous section is to find the lower and upper bounds on the state vector x that satisfy the inequalities (12a)-(12g).To achieve this, a total of 2(n l + n u ) LPs must be solved where each problem will derive either the lower or upper bound of an individual state variable, indicated by x * i .This procedure is Vrachimis et al., August 9, 2017 described in Algorithm 1.
Algorithm 1 Solution of LISE using LP begin Minimize x i by solving problem LPmin 3: Maximize x i by solving problem LPmax 5: x u i = x * i 6: end for 7: x = x l , x u return x Step 5: Iterative solution of the nonlinear interval system of equations In Algorithm 1 the linearized version of the original problem in ( 9) is solved.This is an outer interval solution to the nonlinear problem, which guarantees to include the interval hull solution.To find the smallest possible interval x = x l , x u that satisfies ( 9), an iterative method is used.At each iteration m, the range xbnd in which the optimization algorithm searches for an optimal solution becomes smaller and is redefined as where x(m+1) are the bounds calculated for the state vector x at iteration m.The iterations stop when the bounds on the state vector remain relatively unchanged, i.e. the change is smaller than a small number .The relative change in bounds is computed as follows: The overall procedure for the solution of the NISE is described in Algorithm 2. The resulting bounds are an OI solution of these equations that approximate closely the IH solution.

CASE STUDIES
In this section, the IHISE algorithm is applied on different WDN in order to demonstrate the calculated bounds and evaluate the performance of the algorithm.An illustrative example is given in the "Illustrative example" section, where the bounds on different states Algorithm 2 Iterative solution of NISE begin 1: Define initial bounds h (0) using physical properties.2: Calculate initial bounds q (0) using the procedure in Step 1. 3: x(0) bnd = q(0) h (0) 4: m = 0 5: while e (m) > do 6: Bounding linearization of (9) for x ∈ x (m) bnd 7: Formulate problems LPmin and LPmax 8: Find x (m+1) using Algorithm 1 9: 10: are shown graphically and compared with bounds obtained by MCS.In the section "Results from benchmark networks" a more extensive analysis of the algorithm is presented, as it is applied on different benchmark networks with varying characteristics.The performance of the algorithm is evaluated by defining appropriate performance metrics and by comparing the IHISE bounds with bounds obtained by MCS.In the simulations we assumed a demand-driven modeling approach with no leakages, which translates into f ext,j (h j ) = 1 and q leak,j (h j ) = 0, j ∈ N u in Problem (9).This modelling approach allows us to evaluate the performance of the algorithm using the established WDN simulation software EPANET (Rossman 2000).

Illustrative example
The benchmark network "Net1" shown in Fig. 1 provided by EPANET, is used to demonstrate the bounds on hydraulic states produced by the IHISE algorithm.The network parameters are shown in Table 1.Realistic water demand patterns, are assigned at each demand node.
The IHISE algorithm is used to generate bounds on water flows in pipes and hydraulic heads at nodes of the network.The measurement uncertainty is defined as ±5% on the given water demands at nodes, which is the typical error given by manufacturers of water flow meters.Modeling uncertainty is also considered and it is defined as ±5% on pipe Hazen-Williams coefficients.The simulation duration is 24 hours, with a discrete time step of one hour.
Additionally, the same bounds are generated using Monte-Carlo Simulations (MCS) of the network in EPANET.The demands are randomly varied at each simulation within a range of ±5% of the given water demands at nodes.The uncertainty on pipe Hazen-Williams coefficients is achieved by analogously varying pipe lengths, as the Hazen-Williams coefficients are linearly depended on this parameter.Uncertainty on pipe roughness coefficients and pipe diameter can also be considered, but the effect on Hazen-Williams coefficient will not be linear.The maximum and minimum value of each state is saved, defining the upper and lower bounds.The number of simulations is set to 30 000.Note that MCS provide an inner approximation of the bounds on each state and how close they are to the true bounds depends on the number of simulations.Given the possible variations of the same network for the given uncertainty, a sufficiently large number of simulations need to be performed in order for the MCS to converge to the true bounds.Nevertheless, especially in small networks, the MCS bounds can be useful to evaluate the IHISE bounds, which are an outer approximation of the true bounds, by: 1) verifying the correctness of the IHISE bounds by checking if the MCS bounds are always a subset of the IHISE bounds and 2) evaluating the conservativeness of the IHISE bounds by measuring their distance from the MCS bounds.Simulation results for selected states which reflect the results for all the states are given in Fig. 2. The IHISE bounds are compared with bounds generated using MCS for each state.
The figure illustrates that the MCS bounds are a subset of the IHISE bounds, while they are also closely approximated.Note that, the true unknown bounds are enclosed between the IHISE and MCS upper and lower bounds.

Results from benchmark networks
To evaluate the ability of IHISE to compute state bounds, five benchmark networks with varying characteristics are used: "Anytown", which was used as the basis for the original "Battle of the Network Optimization Models" (Walski et al. 1987), "Net1", "Net2" and "Net3", which are example networks in EPANET (Rossman 2000), and "ky3" from the Kentucky Infrastructure Authority database of water distribution models (Jolly et al. 2014).
The networks and their characteristics are listed in Table 1.
The networks were carefully selected to have multiple varying characteristics in terms of size, topology and types of elements they contain.Varying the network size, i.e. the number of nodes and links, demonstrates the scalability of the IHISE algorithm by considering the performance in networks with different number of states, i.e. heads at nodes and water flows at links.Varying the number of reservoirs and tanks, as well as the number of pumps in the network, reveals the ability of the algorithm to deal with these components.The topology of the networks is also considered, specifically the complexity that arises in calculating hydraulic states when the networks contain loops.This is quantified by calculating the circuit rank of the network indicated by γ, which is then normalized by the number of links n l of each network.The resulting metric is defined as the Loop Ratio, given by LR = γ/n l , 0 ≤ LR < 1.A value of LR equal to zero means that there are no loops in the network, while LR approaches the value of one in the case of a fully connected graph.
Note that the circuit rank of an undirected graph is defined as the number of independent cycles, or the minimum number of edges that must be removed from the graph to break all its cycles making it into a tree.It is calculated as r = m − n + c, where m is the number of edges in the given graph, n is the number of vertices and c is the number of connected components.The circuit rank is also known as the cyclomatic number and is used to indicate the complexity of a program's source code (McCabe 1976).

Monte-Carlo and IHISE Simulations
For each network, a random demand scenario is assigned which produces a feasible solution in EPANET, i.e. there are no negative pressures.Similar to the "Illustrative example" section, the demand and modeling uncertainty is consider equal to ±5%.The simulation duration for all networks is 24 hours, with a discrete time step of one hour.The total number of simulation steps is defined as N s = 24.
MCS were performed for all networks using EPANET.The varying parameters are the nodal demands and pipe parameters, in the range defined by the assumed uncertainty.At each simulation, the minimum and maximum value of each state for each time step is saved, thus defining the lower and upper bounds on the state.Additionally, a robust estimate of the state under the considered uncertainties is calculated by taking all the simulated scenarios and calculating the mean value of each state at each time step.This is indicated by for flow states and h µ M C,j (k) for head states and will be referred to as the MCS state estimate.
In order to perform an appropriate number of MCS and obtain quantifiable bounds for each network, a stopping criterion for the simulations is imposed, such that the change in all upper and lower flow-state and head-state bounds are less than ∆q(m 3 /h) and ∆h(m) respectively, for at least 5 000 consecutive simulations.The flows and heads in each network may belong in different value ranges, so ∆q and ∆h are calculated as a percentage of the absolute mean value of the MCS state estimates.The absolute mean value of flow-states and head-states respectively is given by: The defined accuracies are then calculated as ∆q = 1%(µ q ) and ∆h = 1%(µ h ).Using this approach, it is assumed that the bound accuracy is given by ∆q and ∆h.Note that while this approach provides a degree of confidence for the accuracy of MCS bounds, it is not guaranteed that the deviation of these bounds from the actual bounds cannot be larger.
Using the IHISE algorithm, bounds for the state of all networks are computed, using the assumed uncertainty.As a technical note, due to the fact that the IHISE algorithm was designed from scratch without the use of other hydraulic solvers, the networks had to satisfy specific conditions in order for the current version of the algorithm to work and be compared with the results from EPANET.The tank levels have to be measurable, so any tanks in the network model are replaced with variable head reservoirs.Control rules that open and close pumps depending on tank levels were removed from the model.Additionally, the head-loss formula should be set to Hazen-Williams.These limitations will be removed in future versions of the algorithm.

Evaluation of bounds mean value
An estimated value of the each state at each time step k is derived using the mean value of the IHISE bounds, which we will refer to as the IHISE state estimate.For flow-states this is calculated as q µ IH,i (k) = q u IH,i (k) + q l IH,i (k) /2, and for head-states this is calculated networks, the Mean APD was less than 1%, while 99% of the time, the APD was less than 8%.This indicates that, when no statistical characterization of the uncertainties is available, the mean value of the IHISE algorithm bounds can be used as a robust estimate of the system state.

Evaluation of bounds
For the evaluation of the IHISE algorithm bounds, we compare the lower bounds x l IH = [q l IH h l IH ] and upper bounds x u IH = [q u IH h u IH ] of the IHISE algorithm with the lower bounds First, the validity of the IHISE bounds was checked, i.e. x l IH,i (k) < x l M C,i (k) and x u IH,i (k) > x u M C,i (k) for all states i and all time steps k.The test indicated that there were no bound violations for networks "Net1", "Net2" and "Anytown".For "Net3", bound violations occur for two flow-states, in time steps when the MCS flow estimate is less than 10 −4 (m 3 /h) and the MCS bound width is less than 10 −2 (m 3 /h).For "ky3", bound violations occur in 1.77% of the time, in specific states where the difference in width of the IHISE bounds and MCS bounds is less than 0.4(m 3 /h) for flows and less than 0.005(m) for heads.The violation magnitude for any state and time step is less than 0.5% of the corresponding MCS state estimate.All the observed bound violations occur in cases where the IHISE bounds are very close to the MCS bounds.This can be explained by the fact that the IHISE algorithm uses a completely independent hydraulic solver than EPANET, thus differences in solutions may exist.Despite this fact, the differences in solutions made apparent by the violations are insignificant and the validity test of the IHISE can be considered succesful if these are attributed to modeling uncertainty which has not been taken into account.
Next, the Absolute Deviation (AD) of the two sets of bounds is evaluated separately for flow-states and head-states, as they are measured in different units.The AD for flow-state lower and upper bound is defined as e u q,i (k) = q u IH,i (k) − q u M C,i (k) and e l q,i (k) = q l M C,i (k) − q l IH,i (k) respectively.Similarly, the AD for head-state lower and upper bounds is defined as e l h,j (k) and e u h,j (k) respectively.An illustration of lower and upper bound ADs is shown in Fig. 3.In Table 1 the mean AD for all states and time steps are shown for each network.
The results indicate mean errors for upper and lower bounds that are close to the accuracy of MCS, ∆Q and ∆h.
The area defined by the IHISE bounds is also an important evaluation metric.Since the duration of simulations is the same for all states, evaluating the area is equivalent to evaluating the width of the bounds.The width of the bounds for each time step is defined as the difference between the upper and lower bound for each time step.For IHISE flow-state bounds, for state i and time step k, the width is given by w Similarly, the width of IHISE head-state bounds, for state j and time step k, is defined as w h IH,j (k).The corresponding widths for MCS bounds are denoted by w q M C,i (k) and w h M C,j (k).
An illustration of bound widths is given in Fig. 3 In Table 1 the mean bound widths for all states and time steps are shown for each network.The mean IHISE bound widths is indicated by w q IH for flow-states and w h IH for heads-states, while the MCS bounds are similarly indicated by w q M C and w h M C .As shown in Table 1, the difference between IHISE and MCS mean bound width for flow-states in different networks varies from 0.6(m 3 /h) in "Net2" to 37.5(m 3 /h) in "Net3".Similarly, the mean bound width for head-states varies from 0.04(m) in "Net2" to 2.05(m) in "Net1".
In order for the bounds width to give meaningful insight into the accuracy of the algorithm, they must be normalized relative to the absolute mean value of states for each network, i.e. µ q for flow-states and µ h for head-states.Using this normalization, the bound width can be viewed as a percentage of each state's uncertainty.The percent state uncertainty (PSU) is calculated using the bound width-to-mean ratio as follows: where alg = {IH, M C} depending on the algorithm used, and s = {q, h} depending on the type of state.This will allow the comparison between the calculated state uncertainty and the uncertainty on the network inputs, i.e. the demand and parameter uncertainty.It is recalled that the uncertainty on demands and parameters is defined as a percentage of their estimated value, which was set at ±5% in these simulations.
The average PSU for each network, indicated by η s alg : alg = {IH, M C}, s = {q, h}, is given in Table 1.For flow states, the PSU is, for both methods, close to the ±5% input uncertainty which will be used as a reference point.Typically MCS have slightly less uncertainty and IHISE slightly more, with the exclusion of the looped network "Anytown" where both methods have more uncertainty, and the small network "Net1", where both methods have more.For head-states, the results are much different, as both methods produce much less state uncertainty than the reference point, except in the case of "Net1" where the uncertainty is near ±5%.The IHISE average PSU is at the worst case 2 times larger than the MCS average PSU.The worst cases present at the large network "ky3" but also in the highly looped network "Anytown".The advantage of this methodology over MCS is that the calculated bounds guarantee the inclusion of the true system state, while the iterative nature of the algorithm makes these bounds as tight as possible.An extension of this work is to use the generated bounds for fault diagnosis methods that detect and localize leakages in the network.The proposed methodology can be naturally used with model based fault-diagnosis and robust control methodologies, because many of these rely on the availability of bounding state estimates which are calculated by some knowledge of the system uncertainties.In the case of fault-diagnosis, the bounding state estimates are used to create thresholds, which when violated is an indication of a fault (Puig 2010).Additionally, the bounds on hydraulic states of the network can be used to generate bounds on water quality states, since the dynamics of hydraulic and quality states of a water network are interconnected.
A limitation of this methodology is that it does not model elements whose head function is depended on pressure, such as pressure reduction valves.This is something that will be considered in future work.Other elements that are used in WDN and are not modeled in this work, are pressure control valves, flow control valves etc. Future work will model a variety of additional components to be used with this methodology, and an interval hydraulic state estimation toolkit will be released.Additionally, more extensive simulations on how this methodology deals with pressure-driven demands and pressure-dependent leakages will be provided.

ACKNOWLEDGMENTS
This research is partially funded by the European Union Horizon 2020 programme under grant agreement no.739551 (KIOS CoE), and by the Interreg V-A Greece-Cyprus 2014-2020 programme, co-financed by the European Union (ERDF) and National Funds of Greece and Cyprus, under the project SmartWater2020.SUPPLEMENTAL DATAA detailed description of the derivation of initial bounds on flow states, as described inStep 1 of the IHISE algorithm in Section 3, is provided in Appendix S1.Appendix S2 contains a detailed description of the bounding linearization procedure of Section 3, an illustrative example showing the convergence over successive iterations and illustrative examples of the application of this procedure on nonlinear functions of different hydraulic elements, pressure dependent demand functions and leakage functions.Appendices S3-S5 provide additional simulation data from the Case Study Section.Appendices S1-S5 are available online in the ASCE Library (www.ascelibrary.org).

Fig. 2 .
Fig. 2. Comparison of selected pipes water flow bounds (above) and selected nodes hydraulic head bounds (below), generated by Monte-Carlo Simulations (blue solid area) and the IHISE algorithm (red dashed lines).
(Eliades et al. 2015)2007;Kolev 2004a there are several solutions proposed in the literature that approximate the IH, either with an Outer Interval (OI) solution, which is any interval vector enclosing the IH solution, or with an INner Interval (INI) solution, which is any interval vector that is a subset of the IH solution.Most approaches deal with the problem of finding an OI solution, while the INI solution is used as a measure of the overestimation of the solution(Neumaier and Pownuk 2007;Kolev 2004a).An example of a method for finding the INI solution to a NISE are Monte-Carlo Simulations, where solutions are calculated by randomly generating and evaluating a large number of non-interval equations with parameters within the defined intervals(Eliades et al. 2015).The set of solutions is always a subset of the IH solution.
Finding the IH solution to general NISE is a very challenging problem; even for the general Linear Interval System of Equations (LISE), find-ing the IH k) /2.The IHISE state estimates are compared with the MCS state estimates, by calculating the Absolute Percentage Deviation (APD) of the IHISE state estimates to the MCS state estimates.Note that data from time steps where the MCS state estimate is close to zero are excluded from the evaluation, as they produced large percentages that are not representative of the results.The comparison shows that, for all

Table 1 .
In Appendix S5 of Supplemental Data, an extrapolation of the simulation times for the IHISE algorithm compared to the simulation times of MCS based on the five networks of this case study is given.The estimated simulation time of the IHISE algorithm is always less than the MCS time with the defined accuracy, while the time difference becomes larger for larger networks.The simulation time also depends on the complexity of the network, as it is evident in Table1from the simulation time of the looped network "Anytown".The IHISE algorithm needs more iterations to converge to a solution for this network compared to the "Net2" which has similar number of states but is less looped.Similarly, more MCS are needed for the looped network "Anytown" than "Net2" for them to converge to a defined bound accuracy.
CONCLUSIONSIn this work the problem of estimating bounds on WDN hydraulic states is addressed.A new methodology is proposed that generates interval state estimates.The proposed Iterative Hydraulic Interval State Estimation (IHISE) algorithm generates bounds on hydraulic states of the network, by taking into account the water demand uncertainty and modeling uncertainty in the form of uncertain pipe parameters.The uncertainties are modeled as intervals.The results show that the proposed methodology is able to generate tight bounds on hydraulic states and can be used in place of randomized methods such as Monte-Carlo Simulations (MCS).

TABLE 1 .
Results of the IHISE algorithm on benchmark networks.