Optimal Sensor Placement in Hydraulic Conduit Networks: A State-Space Approach

: Conduit bursts or leakages present an ongoing problem for hydraulic fluid transport grids, such as oil or water conduit networks. Better monitoring allows for easier identification of burst sites and faster response strategies but heavily relies on suffici ent insight in the network’s dynamics, obtained from real-time flow and pressure sensor data. This paper presents a linearized state-space model of hydraulic networks suited for optimal sensor placement. Observability Gramians are used to identify the optimal sensor configuration by maximizing the output energy of network states. This approach does not rely on model simulation of hydraulic burst scenarios or on burst sensitivity matrices, but, instead, it determines optimal sensor placement solely from the model structure, taking into account the pressure dynamics and hydraulics of the network. For a good understanding of the method, it is illustrated by two small water distribution networks. The results show that the best sensor locations for these networks can be accurately determined and explained. A third example is added to demonstrate our method to a more realistic case.


Introduction
Hydraulic models are an essential tool for ensuring safe, reliable, affordable, and continuous delivery of fluids, such as water or oil, to end-users [1]. Due to the physical size and complexity of most conduit networks, the actual operational conditions of these grids are hard to monitor [2]. Besides serving as a digital twin of the network, model simulations of the real system can be used to predict and forecast, in real-time, network flows and pressures under varying hydraulic scenarios, valve configurations, or conduit leakages [3][4][5][6]. In addition, distribution system modelling can help optimize network design, sensor and actuator placement, or facilitate operation of the network through testing of control strategies of pumping and valve configurations [2][3][4][6][7][8][9][10][11]. Models, therefore, act as an active tool for network management and leakage control, instead of the classical passive approach of solely reacting when a defunct asset is detected.
In order to successfully deploy these models, accurate, and up-to-date, insight into the network is required in the form of real-time measurements from sensors placed throughout the network. Sensor placement, operation, and maintenance is costly, meaning there is a tradeoff between network information gain and sensor costs. A wireless sensor network of as few as possible flow and pressure sensors, at key positions within the network, is of vital importance for optimal network management. Therefore, optimal sensor placement poses an ongoing challenge in hydraulic conduit networks.
Various studies have been conducted with the goal of maximizing the diagnostic performance of a system under budgetary constraints by means of applying optimal sensor placement. Current methods often consider optimal sensor placement with regards to leak detectability by simulating hydraulic scenarios with leakages [2,[11][12][13][14]. These studies are mostly based on simulating leakages at various locations within a hydraulic network. The change in each state (pressure/flow), as a consequence of these leakages, is captured in a binarized sensitivity matrix, represented by, e.g., a Jacobian or forward finite differences matrix. The potential sensor locations are then ranked based on the number of burst locations for which significant sensitivities are found [15,16]. This approach can be expanded in several directions, for example, by using artificial bursts achieved by opening fire hydrants in various areas of the network instead of model simulated bursts [9], by considering demand uncertainty [7], by not binarizing the sensitivity matrix [10,17], by using optimization techniques instead of iterative techniques to determine the optimal sensor location [14], or in combination with a leakage localization algorithm [11]. In addition to these studies on hydraulic networks, observability and sensor placement studies on electrical networks have also been performed [18][19][20].
Although practical and efficient, the performance of existing sensor placement techniques, based on virtual leakage simulation, is highly dependent on the accuracy of the hydraulic model. The estimated sensor placement is very sensitive to uncertainty in demand estimates, model parameters, measurement noise, and asset properties, and has to consider different leakage locations and fluid loss rates, in order to provide accurate sensor placement suggestions [21]. Considering all these factors does, not only, result in a high dimensional problem and extensive simulations of hydraulic scenarios but also in high cumulative uncertainties, which exponentially worsens when considering simultaneous placement of multiple sensors. Recent publications, on the topic of tackling this high dimensional and highly uncertain problem, all suggest to invest more research into these uncertain factors and focus on development of smart optimization strategies to reduce the computational load [6,7,11,17,22]. Additionally, an extra source of uncertainty in existing theories is that pressure changes are assumed to take place instantaneously, which, especially, for larger networks is too rough an assumption [11,23].
The objective of this study is to investigate the observability of conduit networks and optimal sensor placement designs, only considering the structure of the hydraulic network model, and without a dependence on dynamic network simulations. By not relying on simulation of hydraulic scenarios, no computationally expensive high dimensional optimization is required. In this study, a linearized hydraulic network model is presented in state-space form, with, as model states, the pressures in all network junctions and the flows through all network pipes. Model outputs are the model states corresponding to network junctions and/or conduits where a pressure or flow sensor is installed. Based on a likely hydraulic scenario, with corresponding stationary network flows calculated with an EPANET model, the original non-linear hydraulic network model was linearized. This linearization step enables conventional observability analysis [24,25] and optimal sensor placement based on maximizing the output energy of observability Gramians [26]. Current research often focuses on smart optimization to reduce the computational load associated with simultaneous placement of multiple sensors [6,7,11,17,22]. This study, however, aims to present an alternative sensor placement framework that starts with state-space modelling. The advantages of using a state-space model, for sensor placement in hydraulic conduit networks, are demonstrated by two illustrative examples and one more realistic example with corresponding EPANET models [27].
We also show how the suitability of each conduit and junction, in the model for sensor placement, can be mapped on a graph of the network. This visualization allows for easy identification of optimal regions in the network for sensor placement. This visual information can be used to combine observability function-based optimal sensor placement with other network-specific knowledge regarding sensor placement. The statespace methodology in this paper has been developed using open source software and is equipped with the capacity to transform EPANET model files into state-space models using a Python 3 algorithm.

Materials and Methods
In order to perform optimal sensor placement in a hydraulic network, based on systems theory, in this study, the network characteristics and dynamics are presented in state-space form. Given the network characteristics, such as conduit lengths, diameters, and roughness, flows and pressures within the hydraulic system can be modelled using continuity and momentum equations for unsteady, nonuniform flow of a slightly compressible fluid in slightly elastic conduits. For each conduit, these assumptions lead to the following set of hyperbolic partial differential equations [1,28]: Here, is the pressure in , is the flow velocity in , is the mass density of the transported fluid in 3 , is the elastic wave velocity in , is the acceleration due to gravity in 2 , is the angle the conduit makes with the horizontal, with the angle taken positive if the conduits slopes upwards in the flow direction, is the Darcy-Weisbach friction factor (dimensionless), and is the diameter of the inside of the conduit in . The distinction between the magnitude of flow velocity | | and directional flow velocity is made to allow for flow in both directions through a conduit. For a thorough observability analysis of systems described by hyperbolic partial differential equations, we refer to [29].
The slope term sin( ) is relatively small for most applications and may be neglected. Even if the slope is taken into account, the term will be interpreted as a disturbance and will therefore not influence optimal sensor placement, which solely relies on flow and pressure dynamics, as well as sensor configurations.
Also, in many applications, the convective acceleration terms ( ⁄ ) and ( ⁄ ) are small compared to the other terms and may, therefore, be neglected [1]. However, in this study, we assume slightly compressible fluid in slightly elastic conduits and thus changes in pressure and flowrate with distance ( ) are not zero.
Furthermore, with = ( + 0 ) and = , Equation (1) can be expressed in terms of the piezometric head = − 0 above a specified level 0 , and volumetric flow rate with conduit's cross-sectional area = 1 4 2 [1,30]. In this state transformation both and are assumed to be constant. The variation of and is still indirectly taken into account by using a finite elastic wave velocity . The elastic wave velocity is a function of various properties of the transported fluid as well as the conduit, but is assumed constant within a pipe, since the changes within a single conduit are assumed small [31]. For water transport without air bubbles through PVC conduits, the wave velocity is estimated as = 1200 [1].
Equation ( , where the Darcy-Weisbach friction factor is a function of the Reynold's number. When solely considering water transport and assuming constant temperature and viscosity, such as is the case in water distribution networks, the friction losses are not dependent on the Reynold's number according to the empirical Hazen-Williams equation [1]. Expressed in volumetric flow rate and piezometric head, the Darcy-Weisbach friction related flow loss We will apply this system of hyperbolic partial differential equation to a conduit network with junctions and conduits, connecting junctions and . Assuming ( ) along conduit may be approximated by ( ∆ / ∆ / ), where is the length of the conduit in the flow direction, we arrive for the head in junction = 1,2, … , at the following equation [32]: The difference between the flow in a conduit at conduit start and end is a consequence of the slight compressibility of the fluid and the slight elasticity of the conduit. For fluid transport, the difference between flow at beginning and end of a conduit is usually very small. Commonly, large pressure changes, as a result of high elastic wave velocity , are assumed immediate. Hence, both the low flow variation within a conduit and the rapid pressure changes motivate the assumption that pressure changes are instantaneous, implying that the head in each junction is always in steady state, thus = 0 and thus ( ) = ( ) ≡ [1]. Consequently, under these assumptions, the dynamics of the system would be solely governed by the momentum Equation (2) Although this equation is very suitable for calculation of hydraulic scenarios and thus, for performing optimal sensor placement through the use of burst simulations, significant and measurable pressure transients do occur [31]. Since modern sensors can operate under sampling frequencies higher than once per second, the damping oscillations, as a consequence of pressure transients and friction in the conduits, can be detected. Since large pressure transients, also referred to as water hammers, can cause conduit wear and bursts, it is important to be able to identify where, how often, and to what extent these transients occur in order to identify their causes and adopt a mitigation strategy. For optimal sensor placement, based on observability analysis, we take into account these transients by assuming  0, and assuming a linear relationship between change in flow rate and flow rate throughout each conduit in the flow direction → : The relative flow gradient in −1 is small, since the compressibility of fluids and elasticity of conduits are very small in most hydraulic conduit networks. As ε is unknown, in this study, we assume it is unknown-but-bounded. Thus, ε is defined on an interval that represents the uncertainty in the values of the compressibility and elasticity. For details about the relative flow gradient, see Appendix A. This assumption (5) yields a system consisting of a linear continuity and a non-linear momentum equation: Notice that Equation (6) is nonlinear with regards to volumetric flow . However, for small perturbations from a specific hydraulic scenario, as a result of steady state computations in EPANET with steady state pressures ̅ and flows ̅ , the model can be linearized around that hydraulic scenario. Thus, we assume | | 0.852 Consequently, the dynamics of the system are then expressed in terms of the 'resistance' Simulation of a steady-state hydraulic scenario is, thus, required in order to obtain estimates for the linearization points ̅ . Although EPANET hydraulic simulations do not include pressure dynamics, these dynamics are still taken into account in the statespace model Equation (7).
Notice that Equation (3) and thus also Equation (7) describe the pressure change as a result of gradients in the flow rates. In steady state, for incompressible fluid and nonelastic networks, the continuity equation at each node is given by = 0. Small deviations in the flow rates through a node, as a result of changing boundary conditions, may lead to small increases or decreases of the pressure in the node, which are also covered by Equation (5). Consequently, as a result of our approximations, for large changes in the hydraulic scenario, new steady states need to be calculated, using, e.g., the EPANET model (Rossman, 2000).
To put these equations in matrix-vector form, we introduce the state vector : = [ 1 , 2 , … , , 1 , 2 , … , ] ∈ ℝ with = + , where the first elements contain the heads , and the remaining elements the flows . We further introduce the output vector = [ 1 , 2 , … , ] ∈ ℝ with = + , where the first elements contain the heads of those junctions equipped with a pressure sensor and the remaining elements contain the flows of those conduits equipped with a flow sensor. In what follows, we assume a pressure sensor is always placed in a junction and a flow sensor halfway on a conduit. Equation (7) can thus be represented in the following form: where ( ) contains the boundary conditions and in what follows = . The dynamics of the system are determined by the system specific parameters , , and , which make up the elements of the  matrix . The positions of the sensors are specified via the  matrix . For applications of optimal sensor placement, only the system dynamics (matrix ) and the sensor locations (matrix ) are required. For model simulations, however, the  input matrix would also be required and would contain inputs such as height differences between junctions, minor losses (valves, pumps), storage junctions (tanks), and the set values of flow or pressure at water sources and sinks (demand or reservoir junctions), and thus, at the boundaries of the system. Thus, for the intended goal of optimal sensor placement, based on state-space methodology, matrices and do not need to be specified and no temporal discretization of the system, Equations (7) and (8), is required.
The output vector is dependent on the junctions with head sensors and the conduits with flow sensors, resulting in a binary pseudo-diagonal output matrix with = + , where those elements of are 1 if a sensor is present at that junction or in that conduit. For any chosen sensor configuration and accompanying output vector and output matrix , the network is observable if the  observability matrix has full rank: A system that is observable allows a full reconstruction of the states over time from given input-output data. However, for large networks, the observability matrix may be ill-conditioned, which would lead to an observability analysis that does not lead to accurate conclusions [33]. If the eigenvalues of matrix all have negative real parts, the system is called (asymptotically) stable. In that case, for each sensor configuration, the  observability Gramian of the network can be calculated, after solving the discrete Lyanupov equation + = − , and is given by: If , a symmetric matrix and unique solution to the discrete Lyapunov equation, is positive definite, that is, has all eigenvalues larger than zero, then the system defined by and is observable. For linear, time-invariant systems, such as given by Equations (7) and (8), the sensitivity of the output with respect to the initial state (0) is given by [33]. Therefore, the observability Gramian from Equation (10) can be interpreted as a Fisher Information Matrix, and it can thus be understood as a measure of information content. Its inverse, apart from a scaling factor, represents the uncertainty in the estimates of the states [34]. In the following, a norm of the observability Gramian will be used as a measure for network observability. In this study, the smallest eigenvalue of is chosen as norm, instead of a summarizing functional based on "optimality criteria" from optimal experiment design [35,36]. Georges showed that the eigenvalue-optimality criterion can be used to determine which system configuration, defined by , as a result of the choice of matrix , maximizes the observability [26]. This is achieved by quantifying the information content or "output energy" ℰ( ) associated with each different sensor configuration, based on the real-valued non-negative eigenvalues of the corresponding observability Gramian : The smallest eigenvalue corresponds to a combination of states that are least observable. Choosing a sensor configuration that maximizes this minimum eigenvalue ensures maximum observability of this combination of least observable network states, thereby realizing the most meaningful increase in network observability. The sensor configuration that maximizes the output energy ℰ( ) is the optimal sensor configuration that maximizes the network's observability: Although the exact magnitude of the observability index, in this case the smallest eigenvalue, of a specific sensor configuration might not be preserved after our approximations, eigenvalue-optimality still allows for comparison of the observability index of different sensor configurations.

Results and Discussion
In order to illustrate the power of the state-space representation of hydraulic conduit networks introduced in Section 2, we will apply the proposed method, for optimal sensor placement, to two hydraulic models of water distribution networks. Since the aim is to show the essence of the method, we restrict ourselves to small networks. In what follows, we assume a constant relative flow gradient = 10 −3 −1 . See Appendix B for an analysis of the effect of the value of on output energy and optimal sensor placement.

Example 1: Triangular Network
The small triangular network we study here is sketched in Figure 1a, and its properties are specified in Table 1. The network consists of three junctions = 1,2,3 that are connected in a loop via conduits = 12,23,13 (Table 2). An additional conduit = 41 connects a reservoir = 4 with constant head 4 0 = 243.84 to node = 1 . The outgoing reservoir flow is assumed to be known, either inferred from measured reservoir volume or directly measured with a flow sensor on conduit = 41. The eigenvalue decomposition of the state matrix of the triangular network is detailed in Appendix C, showing that the system is asymptotically stable. The question, however, is: where could one extra sensor be best positioned?    Since the triangular network is small, the corresponding observability matrix is wellconditioned. Eigenvalue decomposition of the observability Gramian reveals that the three smallest eigenvalues are significantly smaller than the others (Figure 1a). Especially regarding the two smallest eigenvalues, the corresponding weights of the states of node 2 and node 3 are significantly higher than the weights of the other states in the eigenvectors associated with these two smallest eigenvalues. This indicates that the heads in nodes 2 and 3 are significantly less observable compared to the other states, and placing a sensor in either of these nodes will greatly improve the observability of the least observable part of the network. Singular value decomposition of the observability matrix will give the same result, but is less prone to ill-conditioning for large networks and thus presents a more robust indicator for optimal sensor placement.
For observability-based sensor placement, six options, and thus six different matrix, were considered: a head sensor in one of the nodes or a flow sensor in one of the conduits other than conduit 41, since the flow in conduit 41 is already metered. Notice from Equation (10) that the observability Gramian is defined in terms of an inner product. In order to best visualize the output energy differences between various sensor placements, a square root color scale was used to put emphasis on the comparison between the output energies (Equation (11)) of the different sensor placements ( Figure 2). As can be seen from Figure 2, sensor placement in junction 2 maximizes the output energy (smallest eigenvalue), closely followed by junction 3. As discussed above, this is in line with expectations, since nodes 2 and 3 were the states responsible for the smallest eigenvalues of the observability Gramian.  (11)), corresponding with sensor placement in that specific junction or conduit.

Example 2: Net1 Case Study
In order to perform observability-based sensor placement, based on linear state-space models, estimates of the flow ̅ through the network are required, using steady state hydraulic simulation of the system. However, these flows can differ significantly between scenarios. Therefore, an additional investigation was performed to determine the effect of hydraulic scenarios on resulting optimal sensor location. Optimal placement of one additional sensor for the EPANET chlorine decay model named 'Net1′ was also considered (Figure 3a) [27]. Net1 is a network with one reservoir, tank, and pump, where the flow from/to the reservoir and the tank are assumed measurable, either directly or indirectly, from monitoring the reservoir and tank volumes. Depending on the time of day and the tank water volume, reservoir 9 is decoupled from the network, and tank 2 will act as a water source instead of a sink (Figure 3b). Scenarios with and without the reservoir will result in different optimal sensor locations. Therefore, placement of one additional sensor in Net1 was investigated with and without reservoir, at 08:00 and 20:00, respectively. Analysis at intermediate time instants, and thus, for different supplies and demands with corresponding steady state values of and , did show different values of the output energy. However, this did not lead to changes in the optimal sensor location. Placement of one sensor was considered, in addition to the existing flow sensors, at conduit 110 and 9 (if the valve on conduit 9 is open). We found, for the case at 08:00, that a head sensor in junction 31 is optimal regarding network observability, as seen from the maximum output energy of this sensor configuration compared to alternative sensor placements ( Figure 3). Depending on the valve configuration in the network, junction 32 could also be considered for sensor placement (Figure 3b). However, junction 31 is found to be optimal for both network configurations, whereas junction 32 is not significantly more suited for sensor placement compared to 31, and it is significantly less suitable for the network configuration where the reservoir is connected to the network (Figure 3a). In both cases, the optimal sensor location is at the south end (bottom) of Net1. This is to be expected, since the original Net1 network only contains sensors in the north (top) of the network, so an additional sensor in the south enables network-wide insight. The fact that our placement procedure leads to optimal positions that are very close to each other, although the flow conditions are rather different, indicates that the linearization step does not significantly impact sensor placement performance. Since the valve configuration of the network for other time instances is similar to the configurations at 08:00 or 20:00, with only slight differences in network pressures and flows, optimal sensor placement for other time instances yields the same optimal sensor placement results.
If pressure changes are assumed instantaneous, and thus = 0, only Equation (4) would remain for analysis of optimal sensor placement. Consequently, only flow sensor placement will be regarded optimal using this approach. In this case, the best choice is to place a flow sensor in the conduit with the largest resistance. Since both pressure and flow dynamics of the hydraulic system are included in the state-space model (Equation (7)), factors such as pressure wave velocity will effect sensor placement and thus result in more robust placement and investigation of pressure sensor placement in addition to flow sensor placement. In a practical sense, for full real-time reconstruction of all states, thus including the effect of water hammer, high speed (milliseconds-seconds) sampling sensors are needed, which are not commonly used. However, high speed sampling (in the order of milliseconds to seconds) is not considered a challenge nowadays, and the results presented vote for this strategy.

Hanoi Network
The triangular network was used to illustrate the state-space methodology, and the Net1 example network shows the robustness of state-space sensor placement with regards to changes in hydraulic scenario used for linearization. However, both networks are small theoretical networks. In order to investigate optimal sensor placement in a real network, the Hanoi network was used as a third case study. The Hanoi (Vietnam) network, is a drinking water distribution network with 34 conduits and 31 demand junctions, and it is used as a benchmark for optimal network design application [37,38].
Placement of one additional pressure sensor in the Hanoi network was successfully performed using the state-space methodology ( Figure 4). When assuming the reservoir outflow conduit already contains a flow sensor, placement of a single pressure sensor in junction 25 is deemed optimal for increasing the observability of the network's least observable regions. A sensor on the border of the first and second network loop is, therefore, deemed optimal. Since the boundary of second and third loop is already metered via the flow sensor at the network reservoir, this configuration allows for metering all three loops as thoroughly as possible when placing just a single pressure sensor. This, in turn, will result in greatly improving the observability of each region of the network. In practice, this means that placement of one additional sensor will greatly benefit reconstruction of all network states (flows and pressures) and, therefore, will greatly supplement all network-wide methods and models, such as leakage detection algorithms or digital twins based on hydraulic models.

Conclusions
Optimal sensor placement is a well-studied topic within the literature on smart water grids. However, the focus of these studies often lies on reducing placement uncertainty, as well as more computationally efficient optimization of the involved calculations. The method presented in this study expands the often used burst detectability-centric sensor placement criterion to an observability-based criterion. Our approach ensures that placement of additional sensors will provide more information about the entire network and will help improve hydraulic models or digital twins of the water distribution process. Using a state-space approach that takes flow, as well as pressure dynamics, into account, and does not rely on dynamic simulations, optimal sensor placement can be performed with limited computational efforts. Results based on three case studies indicate a robust sensor placement performance solely based on network observability. Additionally, the effect of piece-wise linearization of the system, as a result of changing hydraulic scenarios, is shown to not significantly impact sensor placement.  Fujiwara and Khang, 1990) is available via original publication. The triangular model data presented in this study is fully described in Table 1. The Net1 chlorine decay model is aviable from EPANET [27]. Custom Python 3.7 code was developed for this study.
Acknowledgments: This work was performed in the cooperation framework of Wetsus, European Centre of Excellence for Sustainable Water Technology (www.wetsus.nl, accessed on 7 October 2021). Wetsus is co-funded by the Dutch Ministry of Economic Affairs and Ministry of Infrastructure and Environment, the Province of Fryslân, and the Northern Netherlands Provinces. The authors would like to thank the participants of the research theme "Smart Water Grids" for the fruitful discussions and financial support.

Conflicts of Interest:
The authors declare no conflict of interest.

Notation
The following symbols are used in this paper: Pressure Pa Let us analyze the approximation of flow differences in a pipe, ( ) − ( ) =  ( ) with   ℝ, in some more detail. This approximation is derived as follows. Assume the flow at the end points of a pipe is given by: With  the relative flow change per meter. Let the unsteady, nonuniform flow of a slightly compressible fluid in slightly elastic conduits, after linearization of the friction term around ̅ , be described by the hyperbolic partial differential equation: where resistance ̅ = Hence, choosing  = 1/ will give the eigenvalues or poles of system (A4).

Appendix B
As expected, placing a sensor in the triangular network, in addition to the sensor in conduit 41, will always result in a higher output energy, as an additional sensor will increase system observability, independent of assumed value of the flow gradient ( Figure A1). For values of within the interval [10 −6 , 1], the sensor configuration with a pressure sensor at node 2 or 3 maximizes the output energy. Therefore, in the case studies we chose = 10 −3 −1 , a good estimate regarding the application of optimal sensor placement. Figure A1. Effect of assumption of relative flow gradient on the Output Energy.

Appendix C
The eigenvalues and eigenvectors (row-wise) of the state matrix of the triangular network are presented in Table A1. Notice from the eigenvalues that the system is asymptotically stable and will show oscillatory behavior, as expected. The "fastest" characteristic mode is related to the heads in the three nodes (see 5th row of eigenmatrix). Table A1. Eigenvalues and accompanying eigenvectors of system matrix .