State Variable Implications on Hydraulic State Estimation

Hydraulic State Estimation allows to estimate the most likely state of a water supply network considering multiple sources of information and their associated uncertainty. It is set out as an optimization problem, often addressed according to the Weighted Least Squares criterion. It can be formulated differently depending on the selected set of state variables. This choice is not straightforward and leads to different problem dimensions, time complexities and convergence behavior. All possible essential approaches are gathered in this work and two are identified as the most suitable according to the resulting problem dimensions: a head-based approach and a demand-based approach. The particularities and the formulation according to both sets of state variables are discussed, seeking for efficient implementation. Their time complexity and convergence behavior are compared to draw conclusions that help to identify the most suitable approach for real practice applications.


Introduction
Having a reliable knowledge of the hydraulic state of a water network is essential for its management, which seeks to meet user requirements. Traditionally, water network analysis has focused on models that solve the hydraulics of water systems (e.g. Cross 1936;Martin and Peters 1963;Epp and Fowler 1970;Hamam and Brameller 1971;Wood and Charles 1972;Todini and Pilati 1988). These methods are all projections of the same set of equations (Todini and Rossman 2013): the governing physical equations of pipe flow (i.e. mass continuity and energy conservation) and the network topology. Their differences stem from the driving unknowns adopted to solve the problem. Note that the usual inputs required to solve the flow problem in a basic network are head levels at knownhead nodes (reservoirs and tanks), nodal demands and pipe roughness values. Pipe flows and the rest of head levels constitute the problem unknowns (Elhay and Simpson 2011). Either head levels or flows can be considered the driving unknowns, leading to two different approaches to the problem: computing head levels and updating flows or computing flows and updating head levels. The Global Gradient Algorithm -GGA method-(e.g. Todini and Pilati 1988;Rossman 2000) and the Loop method (e.g. Epp and Fowler 1970;Alvarruiz et al. 2018) are good representatives of these head-based and flow-based approaches, respectively.
Regardless of the method adopted, the hydraulic model must be able to simulate the complex reality of the water system. In practice, not only traditional inputs (head levels at known-head nodes, nodal demands and roughness values) but any system variable may be measured or estimated based on historical records (i.e. pseudomeasured). Moreover, these measurements and/or estimates are always uncertain (Bragalli et al. 2016), as different sources of information may be subjected to different levels of uncertainty (Ruiz et al. 2022). Therefore, it can be argued that the traditional scheme of solving the flow network has two substantial drawbacks when coping with the complex reality of water systems: (1) the input set is closed, so other network measurements cannot be readily incorporated and (2) inputs do not systematically consider uncertainty in each model run. These two setbacks can be sorted out with simulation methods (Kumar et al. 2008), but they are computationally expensive because they need several model runs to characterize a single time.
Hydraulic State Estimation (HSE) can be considered an alternative to overcome these difficulties. HSE is posed as an optimization problem that respects the subjacent flow governing equations (Carpentier and Cohen 1991). All water system variables (i.e. head levels, nodal demands, pipe roughness values and pipe flows) are considered unknowns in this approach. They must be computed considering the set of measured variables and the set of state variables (Jung and Kim 2018). Measured variables come determined by the measurements and / or pseudomeasurements available at the network (Bargiela and Hainsworth 1989). This enables to consider any measurement/pseudomeasurement and their uncertainty as input. State variables are the group of independent variables required to describe the network hydraulic state, i.e. to compute the rest of the system variables (Sterling and Bargiela 1984). State variables are related to the rest of variables by the governing physical equations of pipe flow and the network topology (Bargiela 1984). Therefore, HSE solves the hydraulics of the water system, but it goes one step further: it provides the value of the state variables that minimizes the difference between computed measured variables and measurements/estimates (Díaz et al. 2017).
As it happened when solving the flow network problem, several approaches can be posed to address the HSE problem (e.g. Powell 1992;Andersen and Powell 2000;Díaz et al. 2017Díaz et al. , 2018Vrachimis et al. 2019;Wang et al. 2021). These approaches are characterized by the objective function criterion and the selected set of state variables. Weighted Least Squares (WLS) is one of the most popular criteria for the objective function in HSE literature due to its proved advantages (Tshehla et al. 2017). Regarding state variables, to the best of the authors knowledge, there is no thorough discussion about the implications of adopting a specific set of state variables. The choice is not trivial from a computational point of view: different sets lead to problems with different dimensions, convergence performance and time complexity. Hence, it is important in order to operationally implement HSE in large real systems.
The objective of this paper is to provide a conceptual framework to assess the computational implications of selecting a specific state variable set for HSE and to identify the most suitable approach in real practice. For this purpose, time complexity, execution time and convergence performance are here evaluated. An efficient computational implementation of the problem is proposed and formulae are developed for any possible measurement at a water system. HSE is formulated all along according to the WLS scheme because it is the most widely adopted criterion in the literature and it presents several advantages that will be later explained. A similar state variable analysis could be undertaken for other objective function criteria if needed, but it is out of the scope of this paper.
The rest of the paper is organized as follows. An overview of HSE under a WLS approach is first presented. Then, all possible essential sets of state variables are considered and two (a head-based and a demand-based approach) are identified as the most efficient, bearing resemblance with the two options that have been mainly adopted to solve the traditional flow network problem. These two sets are discussed in more detail, explaining their implications on the formulation and giving hints for efficient implementation. Their time complexity is generally analyzed and execution times and convergence performance are assessed for several case studies. Finally, conclusions are drawn in order to identify the most suitable approach.

Overview of Hydraulic State Estimation
HSE is an optimization problem that provides the most likely state of a water supply network. This state is understood as the one which minimizes the difference between the computed measured variables and their corresponding measurements or pseudomeasurements. Such variables are calculated from the estimated state variables. Both state and measured variables are related non-linearly through the governing equations of water systems (Díaz et al. 2016b), i.e. mass continuity, energy conservation and network topology. This is expressed with the following model: where ℝ m is the vector of measurements, ℝ n is the vector of state variables, ( ) ∶ ℝ n → ℝ m is the function of non-linear relationships between measured and state variables and ℝ m is the vector of differences between measured variables and their corresponding measurements, which is to be minimized. It is important to note that ( ) is a subspace of ( ) ∶ ℝ n → ℝ v , which is the function of non-linear relationships between state variables and any other water system variable.
Several criteria might be chosen to minimize the model in Eq. (1), leading to different objective functions. State estimation was first adapted to water supply systems by Coulbeck (1977). He used a WLS approach for the objective function, although measurement weights were arbitrary. This criterion was chosen because the unweighted Least Squares yields the minimum variance unbiased solution when measurement uncertainties are solely Gaussian (Tshehla et al. 2017), but water system uncertainties are not and they register significant outliers (Bargiela 1984). Measurement weights are used to deal with high uncertainties while maintaining the virtues of the unweighted Least Squares. Bargiela (1984) defined measurement weights as the inverse of measurement uncertainty, establishing a clear relation between those two. He also compared WLS method with (1) = − ( ) the Weighted Least Absolute Values criterion. His studies drew that both approaches were comparable in terms of numerical stability and convergence, but WLS was a better option for large, hardly measured and noisy systems, such as water supply networks. These advantages explain why WLS has been widely adopted in the literature (and so in this work) for HSE. The objective function according to WLS criterion is posed as (Carpentier and Cohen 1991): where J( ) is the objective function to minimize, ℝ m x m is the variance-covariance matrix of measurements, which acts as the weight matrix, and ( ) ∶ ℝ n → ℝ c is the function of equality constraints. It is important to note that (1) is typically assumed diagonal, namely measurements are independent, and (2) ( ) is also a subspace of ( ) . To specify the constraints, it is necessary to select the set of state variables.
The optimal solution of Eq.
(2) and Eq. (3) is ̂ ℝ n . As this solution minimizes the residuals between measurements and variable estimates, it represents the most likely hydraulic state of the water system under study. HSE allows to get distributed information from scattered measurements and/or estimates, as shown in Fig. 1. In this figure, several disjointed measurements and pseudomeasurements (i.e. inputs) are available, each one with a given uncertainty. From these inputs, distributed hydraulic state information is estimated, filling unmeasured gaps. Since inputs are uncertain, the most likely hydraulic state does not necessarily coincide with every input.
Several numerical methods are available to find ̂ . They consist on iteratively moving from one state to another through the objective function domain. This displacement might be accomplished in different ways and many of the methods achieve it by computing a decrease direction ( ℝ n ) per iteration (Bazaraa et al. 2006), as shown in Fig. 2. Newton's method is widely used in literature since J( ) is twice differentiable and ( ) derivatives can be computed analytically: where i is a counter for iterations, is the vector of total variation of each state variable per iteration, ∇ is the gradient operator with respect to the state variables and ∇ 2 is the Hessian operator with respect to the state variables. It is a gradient method, i.e. the decrease direction is wisely computed using information about the direction of steepest Fig. 1 HSE effect in a water system: isolated data vs distributed information descent. An advantage of Newton's method is that its convergence is quadratic, unlike gradient methods that only use first derivative information, whose convergence rate is linear (Bazaraa et al. 2006).
As shown in Eq.
(2) and Eq. (3), the problem is actually constrained. To take into account the constraints, Newton's method is not only applied to J( ) but also to the Karush-Kuhn-Tucker first order optimality conditions. Thus, i+1 computation takes the form of (Gómez-Quiles et al. 2013): where ℝ m x n and ℝ c x n are the Jacobian matrices of ( ) and ( ) with respect to the state variables, respectively, and ℝ c is the dual variable vector. and are sub-matrices of ℝ v x n , which is the Jacobian matrix of ( ) with respect to the state variables. The necessary and sufficient condition to be able to solve the problem is that has full rank, i.e. the water system is observable (Díaz et al. 2015). This can be achieved if there are enough measurements, estimations and constraints to estimate every state variable. Redundant measurements are recommended because that leads to better system matrix conditioning (e.g. Bargiela 1984).
In each iteration, once is known, the rest of variables of the corresponding hydraulic state are computed through ( ) . This last step is equivalent to solving the flow network according to the governing equations of pipe flow. HSE is not independent from solving the flow network: it embeds and broadens the plain resolution of the physics of the problem. The resulting linear system in Eq. (5) is solved iteratively starting from some initial state variable values ( 0 ) until the stopping criterion is met. This typically corresponds to || i+1 || being lower than a chosen tolerance. To ensure constraint accomplishment, || ( i+1 )|| being lower than a chosen tolerance might also be considered. The algorithm is summarized in Algorithm 1.

Fig. 2 Iterative displacement in the HSE objective function domain
It is important to note that, in Eq. (5), the Hessian of J( ) is approximated. It is assumed that the second derivative of ( ) is null. Caro et al. (2011) state that this approximation is fair enough and holds whether is diagonal or not. However, in order to correctly apply Newton's method, it is important to assume both the Hessian approximation and diagonal structure. This assures that the Hessian is symmetric positive definite. If the Hessian is not positive definite, it cannot be guaranteed that the computed decrease direction is an actual descent direction to the optimum in J( ) domain (Castillo et al. 2002).

Possible Essential Sets of State Variables
Until now, the presented formulae are quite general and could be easily applied to any state estimation problem (i.e. not necessarily a water system). The key point to make the formulation hydraulic-specific is the function ( ) , and its subsequent mathematical objects: ( ) , ( ) , , and . They depend on the selected set of state variables. As stated in Sect. 1, this set is not unique and several configurations are possible. All essential possibilities for a basic network are gathered and compared in Table 1.
ℝ r+j represents head levels, ℝ p stands for pipe flows, ℝ j is nodal demands and ℝ mat is used for pipe material roughness. Regarding dimensions, r, j, p and mat stand for known-head nodes, junctions (rest of nodes), pipes (links) and pipe materials, respectively. can then be split into ℝ r for known-head nodes and ℝ j for junctions. Pipe material roughness (i.e. roughness of each material) and pipe roughness (i.e. roughness of each pipe) are here distinguished to consider that several pipes can be associated with a same roughness value. Both are related through ℝ mat x p , which is a matrix that links pipe material roughness with corresponding pipes. Its components are 1 if a pipe is made of a material and 0 otherwise. ℝ p is the vector of pipe roughness, obtained from and . If mat = p , then = and = . Table 1 shows that different sets of state variables lead to different problem dimensions. It is expected that the higher the problem dimension, the slower the computational time. In water supply networks, j ≤ p and mat ≤ p necessarily (mat is usually much smaller than p) and r ≪ j also holds in practice. It is important to note that, for clarity, the sets presented in Table 1 are the essential possibilities. They can be combined leading to plenty of secondary sets, but this does not provide any advantage in terms of dimension. For example, set has a higher dimension than set and both can be combined resulting in a mixed set with flows and demands. Any nodal demand can be swapped by all flows entering and leaving the node and vice-versa. The dimension of the resulting set is always between the dimensions of and sets. Then, the mixed set has no advantages according to dimensions because set can be used instead and always has a lower dimension. This reasoning can be applied to all combinations.
From the state variable sets in Table 1, two are of special interest: , and (headbased), and , and (demand-based). These two sets are the ones associated with the lowest problem dimension, so they are the most computationally efficient. These two possible sets of state variables for HSE (head-based vs demand-based) bear resemblance with the choice between head levels or flows as driving unknowns when just solving the flow problem.

Head-Based State Estimation
In this case, state variables are , and . With these state variables, ( ) can be defined in three blocks as: where ℝ p is the vector of pipe head-losses, ̂ ∶ ℝ p+mat → ℝ p is the Hazen-Williams head-loss formula and ℝ (r+j) x p is the incidence matrix of the water system topology. components are 0 if a node is not connected to a link, 1 if a node is the final end of a link and -1 if a node is the initial end of a link. is split into ℝ r x p , only considering known-head nodes, and ℝ j x p for junctions. The second block of Eq. (6) fully depends on the selected head-loss formula (Hazen-Williams in this work) and non-linearity stems from it.
It is important to note that, with this state variable set and ( ) defined as in Eq. (6), the rest of system variables (head-losses, flows and demands) are explicit by blocks in ( ) , and these blocks can be sequentially and directly solved. Then, once state variables are available (i.e. computed in each HSE iteration), the rest of variables can be computed straight away. Traditional iterative numerical methods to solve non-linear systems of equations , and r + p + mat , and r + p + j , and r + j + mat , and r + j + mat , and r + j + j , and p + j + mat , and j + j + mat (typically used when solving pipe flow governing equations, e.g. in Loop or GGA methods) are not required to update the hydraulic state through ( ) . This gives computational enhancement to the algorithm. Junctions can be divided into transit nodes ( j t ), with nil demand, and demand nodes ( j d ). In this work, transit nodes are considered fixed constraints that impose null demand. This is the only type of constraint with this set of state variables ( = ; ℝ j t ). According to Eq. (6), these are non-linear constraints with respect to the chosen state variables. Since this division affects demands and not head levels, state variables must still account for all junctions.
According to the state variables and Eq. (6), the structure of is: where is the corresponding dimension identity matrix and ℝ p x p are diagonal matrices coming from the following derivatives: = ∕ and = ∕ . and are obtained from Eq. (7). Adequate matrix rows must be selected according to applicable constraints and available measurements. In case of full measurement redundancy (i.e. all water system variables are measured), all rows are selected either as constraints or measurements. In reality, financial resources are limited and not all possible variables are measured (Savic et al. 2009). Matrices are presented in this work for the full redundancy case for the sake of completeness.

Demand-Based State Estimation
In this case, state variables are , and . This set poses a problem regarding computation. Every derivative ∕ can be obtained through ∕ according to the chain rule. The issue is that the relation between and [ = ( ) ] is not always bijective. To see this, one might think of a looped water supply network. In that kind of network, several sets of are related to the same image if energy is not considered. A necessary and sufficient condition to be able to invert a function is that it must be bijective (Devlin 2004). Then, = ( ) cannot be inverted to = ( ) and ∕ does not exist in this case.
To solve this problem and be able to compute , it is necessary to resort to an equivalent water system model. The function = ( ) is bijective only when the water system is a tree network (in these networks, sets unequivocally and vice-versa without considering energy). Working with the spanning tree of a looped network and considering chord link flows ( ℝ L ) as additional state variables is a valid solution to compute . In this case, ( ) blocks can be redefined as: ℝ T is the vector of tree pipe flows, L is the number of chord pipes or loops (L equals p − j ), T is the number of tree pipes (T equals j), ℝ L x p is the loop matrix, H r 1 is the head level of the known-head node contained in the tree network, ℝ j is a vector whose components are 1 and , ℝ j x T is the sub-matrix of when only tree pipes are considered. For each chord pipe, it is necessary to find a related loop in the tree network. associates loops and pipes. Its components are 0 if a pipe does not belong to a loop, 1 if a pipe belongs to a loop and both pipe and loop directions coincide and -1 if a pipe belongs to a loop but pipe and loop directions do not match. Again, the second block of Eq. (8) fully depends on the selected head-loss formula (Hazen-Williams in this work) and non-linearity comes from it.
Equation (6) could alternatively be used as ( ) with this set of state variables but, in this case, it would be necessary to resort to iterative numerical methods to solve the non-linear system of equations (e.g. as in Loop or GGA methods) in order to update the hydraulic state in each iteration. Redefining ( ) as Eq. (8) for the demand-based approach allows to maintain the advantages of solving the unknowns (flows of tree pipes, head-losses and heads of junctions) explicitly by blocks. Thus, traditional iterative numerical methods to solve non-linear systems of equations are neither necessary with these state variables. This gives computational enhancement to the algorithm.
There are two types of constraints in this demand-based approach. The first type is again null demands at transit nodes, which are here considered as constraints that impose a fix demand ( = ). The second type comes from the need for imposing the physical constraints of the original looped network, which would not be present otherwise in the tree network. Nodal head levels at chord pipe ends must agree with pipe head-losses [ =̂ ( , ) = − , − , ]. Here, , ℝ j x L and , ℝ r x L are the sub-matrices of and when only chord pipes are considered, respectively, ℝ L is the vector of chord pipe head-losses and ℝ mat L is the vector of chord pipe material roughness.
The first type of constraint directly affects some state variables ( ), so it can be imposed by straightly considering those state variables as null. Thus, the amount of state variables to estimate reduces. The resulting state variables are , ℝ j d , and . According to Eq. (8), the second type of constraint is non-linear with respect to the chosen state variables.
According to the state variables and Eq. (8), the resulting structure of is: , related to demand junctions and the other matrices are: where dimensions are: ℝ mat x T are the sub-matrices of , , and when only tree pipes are considered.
ℝ p x p is a diagonal matrix coming from ∕ and ℝ T x T is its sub-matrix for tree pipes. and are obtained from Eq. (9). Adequate matrix rows are selected and combined according to available measurements and energy conservation at chord pipes (second type of constraint). Like before, all rows would be selected in case of full measurement redundancy, which is unlikely to happen in real water supply networks.

Pre and Post-Processing
To perform HSE, pre-processing might be necessary. From Sect. 2.4, it can be guessed that demand-based HSE needs a mandatory topological pre-processing. A spanning tree of the network and the related set of chord pipes and loops must be previously identified. The two most common and broadly used algorithms to find a spanning tree are Kruskal and Prim's algorithms (see Cormen et al. 2001). With respect to loops, a widely used method consists on obtaining a spanning tree of the network and adding loops connecting both ends of chord pipes through the tree. According to Alvarruiz et al. (2015), this method does not suit the typical Loop method to solve pipe flow equations because resulting loops have high overlapping. However, the method behooves HSE because it can benefit from the necessary spanning tree computation and chord links consideration as state variables and overlapping is not a critical issue. Alternative methods might be the ones proposed by Creaco and Franchini (2014) or Alvarruiz et al. (2015). When using head-based HSE, topological pre-processing is not mandatory, although it can be used to find looped and tree zones in complex networks.
Note that, in general, Eq. (9) components are more complex to obtain than those of Eq. (7). Part of this additional computational complexity comes from purely topological matrix operations. These operations can be performed as a pre-processing since their result is not going to change unless network topology changes (e.g. when there are control elements in the water system).
Post-processing can be an excellent add-on to make the most of HSE. The uncertainty of both the estimated state variables and the rest of computed water supply network variables can be assessed. First methods for uncertainty assessment were based on Monte Carlo and Latin Hypercube approaches (e.g. Pasha 2006;Kang 2009). These approaches require = −( , , ) − , many simulations. As one of the main advantages of HSE is that it estimates the most likely state in one simulation, First-Order Second Moment approaches (Caro et al. 2011;Díaz et al. 2016a) may be more convenient to assess HSE uncertainty. Post-processing is here mentioned for the purpose of motivation and completion, but it is out of the scope of this paper.

Results
Two possible sets of state variables have been identified as the most advantageous for HSE based on the problem dimension (i.e. number of state variables). Both of them can be applied to any water supply network, so it is important to assess their performance in order to choose the most appropriate in each situation. Two indicators are used to compare the behavior: time complexity and convergence performance. Time complexity is here expressed in Big-O notation using a uniform cost model, i.e. operations are assumed to have a constant temporal cost. Big-O notation gives an asymptotic upper bound to computational time increment when the input dimension increases (Cormen et al. 2001). Big-O notation is not strictly computational time, so execution times are also analyzed in this work for several case studies. Convergence performance refers to convergence speed (i.e. the number of iterations needed to converge). Convergence is here checked for the same case studies. Bad time results may be counteracted by a good convergence performance and vice-versa.
When assessing Big-O notation, worst-case scenarios and large input dimensions (i.e. close to infinity but computationally manageable) must be assumed to characterize the upper bound. Since Big-O notation describes asymptotic behavior, two approximations are usually taken into account when adding time complexities (Knuth 1997). First, if the addition involves time complexities of the same order, resultant numerical coefficients are neglected. Second, if the addition involves time complexities of different orders, the result is the time complexity with the highest order.

Pre-Processing
Pre-processing is mandatory in demand-based HSE and optional in head-based HSE. Since it is the worst-case scenario for time complexity, it is here assumed that head-based HSE includes an equivalent topological analysis to that of demand-based HSE. As commented in Sect. 2.5, demand-based HSE pre-processing has two stages: topological analysis (tree and loops computation) and topology-related operations of (to speed calculations of matrices in Eq. 9). Pre-processing time complexities are shown in Table 2.
Regarding the spanning tree computation, both Kruskal and Prim's algorithms are well known and their time complexity has been vastly studied. Unless Prim's algorithm is combined with sophisticated heap queue data structures (whose practical advantage is disputed), both algorithms have the same time complexity (see Brodal 1996;Cormen et al. 2001;Brodal et al. 2012). In order to find loops, it is supposed that they are obtained from chord pipes and the previously computed spanning tree, as commented in Sect. 2.5. Regardless of the method, loops are built using shortest path algorithms. For unweighted graphs (such as water systems), a widely used and well known option is the Breadth-First algorithm (Cormen et al. 2001), which is here adopted. Concerning topology-related operations, their time complexity considers taking advantage of the highly sparse structure of the involved topological matrices.
As shown in Table 2, time complexity of demand-based HSE pre-processing is dominated by topology-related operations and it tends to O((j d + j t ) 2 ) , i.e. it only depends on the amount of junctions. By contrast, time complexity of head-based HSE pre-processing is dominated by loop finding and tends to O(L × (r + j d + j t )) . It depends on the number of nodes and loops. While L < j , which is usual in real water systems, time complexity of head-based approach tends to a linear rate and it is smaller than time complexity of demand-based approach, which tends to a quadratic rate. Otherwise, both time complexities are similar if j ≃ L and demand-based time complexity is smaller than head-based time complexity if j < L.

Vector Assembly
This time complexity comprises vectors and assembly and matrix − computation. This time complexity is presented in Table 3 per iteration, i.e. expressions consider only one iteration of the while loop in Algorithm 1. The assembly of vectors and directly depends on the number of state variables and measurements, respectively. It is here assumed all along that measurements have full redundancy, i.e. all water system variables are measured, since it is the worst-case scenario for time complexity. With respect to − , a diagonal matrix structure is assumed, i.e. measurements are considered independent among each other, as discussed in Sect. 2.1.
Both demand-based and head-based HSE tend to the same linear rate, which is dominated by measurement-related operations in case of full-redundancy: O(r + j d + j t + mat + p) . It depends on the total amount of water system variables.

G Computation
Time complexity of computation is determined by those derivatives which (1) imply matrix operations and (2) cannot be completely pre-processed because they are not only topological. Results are again shown per iteration in Table 4. Note that time complexities take Table 2 Pre-processing time complexity

Operation
Head-based Demand-based into account the structure of the involved matrices and which matrices are pre-processed, as well as the operation order. As it will be discussed later, operations according to the headbased approach involve highly sparse matrices, while matrices used in the demand-based approach are highly dense. Time complexity of demand-based computation tends to O((j d + L + mat) × (j d + j t ) 2 ) . It depends on the number of junctions and loops. Time complexity of head-based computation tends to O((j d + j t ) 2 ) , i.e. it only depends on the amount of junctions. The time complexity of head-based approach is quadratic and smaller than that of the demand-based approach. Depending on the relationship among loops, transit nodes and demand nodes at a particular network, demand-based time complexity can be closer to a quadratic or a cubic rate.

Newton's Method
Applying Newton's method involves three major operations: assembling the matrix and right-hand vector of Eq. (5), solving the obtained linear system of equations and updating the state variables with Eq. (4). The time complexities of these operations are gathered per iteration in Table 5, taking into account the structure of the involved matrices. Regarding the head-based approach, both and the matrix of the system of equations have sparse structures. On the contrary, their structures are dense for the demand-based approach.
Under the assumptions in Sect. 2.1, the system matrix of Eq. (5) is symmetric positive definite. Then, plenty of direct methods are available to readily solve the system of equations. Usually, these methods have a cubic time complexity for dense matrices (e.g. see Golub and Van Loan 1996) but they can go quadratic for sparse matrices.
Both demand-based and head-based HSE tend to time complexities which depend on the number of state variables. Head-based time complexity tends to a quadratic rate and it is smaller than demand-based time complexity, which tends to a cubic rate.

Updating Hydraulic State
This time complexity implies computing the three blocks of Eqs. (6) or (8) once state variables have been updated. As mentioned in Sect. 2, Eqs. (6) and (8) can be solved explicitly by blocks. Resultant time complexities are shown in Table 6 per iteration. Table 4 computation time complexity Operation Head-based Demand-based Both demand-based and head-based HSE tend to a linear rate. Until now, numerical coefficients have been discarded when showing time complexity total results, as discussed previously and shown in Table 6. Nevertheless, time complexities of each block and between head-based and demand-based approaches are now very similar, i.e. neglecting numerical coefficients may lead to misleading comparisons. When numerical coefficients are taken into account, head-based time complexity is O(2p + r + j d + j t ) and demandbased time complexity is O(2p + 2(j d + j t )) . It can be seen that the head-based rate is slightly smaller than the demand-based rate. If asymptotic approximations are considered and numerical coefficients are neglected, as in Table 6, both head-based and demand-based time complexities are virtually equal. These approximations are still valid, but the issue is raised here only to highlight that the head-based rate is slightly more advantageous.

Execution Times for Several Case Studies
According to previous results, head-based HSE has a more advantageous time complexity. The benefits of head-based time complexity are significant in pre-processing, computation and Newton's method. These results are now tested comparing the execution times of 6 different case studies: Modena, KL, Sabsevar, Neyshabur, Torbat and Kashmar. They are close to real, extensive water systems that can be used to derive practical conclusions in a broad range of water networks. Important information about the topology of each case study is presented in Table 7. It can be observed that there are less loops than junctions in all seven case studies, as discussed when pre-processing time complexities were explained.
Modena (Bragalli et al. 2012) case study is available at the Centre for Water Systems of the University of Exeter. KL case study is a modified version of the network presented by Kang and Lansey (2012), available at ASCE Task Committee on Research Database for Water Distribution Systems. Modena and KL case studies are well known in water system analysis literature. Sabsevar, Neyshabur, Torbat and Kashmar case studies are generated upon Iranian cities using DynaVIBe-Web application (see Mair et al. 2014a, b) and they are available at its repository. Sabsevar, Neyshabur, Torbat and Kashmar case studies here presented have not been studied in literature before. Layouts of case studies are presented in Fig. 3.  Execution times are shown per iteration in Table 8 for both demand-based and headbased HSE. Measurements are generated from prior solutions from EPANET (Rossman 2000) for each case study. EPANET is a benchmark software in water system analysis and it implements the previously mentioned GGA method (i.e. head-based flow network resolution). In order to resemble reality, only some EPANET's results are used as measurements, i.e. there is not full measurement redundancy, to simulate measurement scarcity. Since it is the usual case in real practice, it is supposed that either measurements or pseudomeasurements of water demands, roughness values, head levels at known-head nodes and flows leaving or entering known-head nodes are available. Only some of the remaining flow and head-level possible measurements are considered to be available. They are randomly taken per case study using a binary random number generator, i.e. generating 0 and 1 with equal probability. Noise is introduced to pipe flow and junction head level measurements to model uncertainty. It is assumed that possible values of noisy measurements follow a normal distribution with EPANET's results as mean values ( ) and given standard deviations ( ). Standard deviations of 1 m and 7.5 are here considered for head levels and pipe material roughness, respectively. A coefficient of variation ( / ) of 1.1 is set for pipe flows and nodal demands. In these tests, noisy values might be within the 95% confidence interval ( ± 2 ). The tests are carried out using MATLAB R2020a version on a computer with 16 GB of RAM and 1 Intel(R) Core(TM) i7-10510U processor at 1.80 GHz and 4 cores.

Solving system of equations
From Table 8, it can be observed that execution times support the obtained time complexities. Execution times are related to problem dimensions and head-based HSE is faster than demand-based HSE. One of the main reasons for the time difference is the distinctive structure of both and the matrix of the system of equations according to the two approaches. These structures are shown in Fig. 4 for Kashmar case study. In this figure, only non-zero matrix elements are colored in grey. It can be appreciated that matrices are sparser under the head-based approach. Moreover, it was stated in Sect. 2.5 that loop overlapping is not a critical issue in HSE. The part of demandbased that is influenced by loop overlapping is highlighted in black in Fig. 4. This portion is small compared with the rest of the matrix.

Convergence Performance
Convergence performance of both head-based and demand-based HSE is checked and compared in the same 6 case studies used before. Convergence performance refers to convergence speed (i.e. the number of iterations needed to converge) in this work. In these tests, each approach is tested under two different assumptions in order to check its robustness against measurement uncertainty: noisy values might be within the 95% confidence interval ( ± 2 ) or within the 99% confidence interval ( ± 3 ).  Modena  268  317  23  49   KL  934  1273  311  339  Sabsevar  9421  9702  6160  281  Neyshabur  9577  9883  6584  306  Torbat  7789  7996  4937  207  Kashmar  4383  4592  2244  209 Tolerances for convergence adapt EPANET's default tolerances (Rossman 2000). The convergence tolerance for the algorithms presented in this work is the norm of being, at least, 3 orders of magnitude lower than the smallest component of and smaller than 10 −6 . Furthermore, the admissible error when imposing the problem constraints is also set. The norm of the vector of errors associated with constraints must be smaller than 10 −4 and, at least, 2 orders of magnitude lower than the smallest value of the variables involved in the constraints. Results are shown in Fig. 5. From Table 7 and Fig. 5, it can be appreciated that, in contrast to time complexity and execution time, convergence performance does not appear to be related with problem dimension. Convergence seems to have a strong relationship with the particular  same iterations. The convergence speed difference does not maintain a constant proportion among case studies. Note that, in the case studies where the demand-based approach converges in fewer iterations, the convergence speed difference is not high enough to counteract the slower execution time per iteration of this approach against the head-based approach. Secondly, demand-based HSE convergence appears to be more robust to high uncertainty. The number of iterations needed to converge in the demandbased approach are the same under both the 2 and the 3 assumptions. In the headbased approach, the iterations needed to converge under the 3 assumption increase with respect to the ones needed under the 2 assumption in 3 out of the 6 case studies.

Conclusions
HSE is a suitable method to estimate the most likely state of a water distribution network. It can overcome some disadvantages of typical flow network resolution models (e.g. GGA or Loop methods) regarding input redundancy and uncertainty without resorting to computationally expensive simulation methods (e.g. Monte Carlo simulations). HSE can be addressed with different optimization criteria, being WLS profusely used in the literature due to its benefits. WLS supports several formulation possibilities depending on the selected set of state variables. Choosing the set of state variables is not straightforward since different sets lead to problems with different dimensions, convergence performance and time complexity. In this work, all possible essential sets of state variables are gathered. Two of them are identified as the most efficient since they are related to the lowest possible dimension of the problem. One of the sets is based on head levels and the other one is based on water demands. The particularities of the so-called head-based and demandbased HSE are presented in this paper together with a computationally efficient formulation. The functioning of both sets is also tackled: their time complexity and convergence performance are addressed. Both obtained time complexities and execution times clearly favor head-based HSE. Head-based HSE has a more advantageous time complexity, especially regarding preprocessing, computation and Newton's method. The obtained execution times per iteration of the tested case studies are much faster in the head-based approach. This is greatly induced by the sparser structure of the matrices involved in the head-based approach. Since time results have been mainly assessed per iteration, the convergence performance (i.e. number of iterations) is also important. Convergence behavior appears to be case specific, although the demand-based approach seems to converge with less iterations, especially for great uncertainties. The lower number of iterations is not able to counteract the worse execution times per iteration.
Finally, it is important to highlight that this paper constitutes an effort to systematically assess the implications of the state variable choice for HSE. According to the results, headbased HSE seems to have a faster performance. It usually needs a few more iterations to converge than demand-based HSE, but its lower per-iteration execution times compensate this difference. Hence, Head-based HSE has been here identified as the most advantageous approach to undertake HSE according to a WLS criterion in real practice.