# First and Second Order Sensitivities of Steady-state Solutions to Water Distribution Systems †

Olivier Piller ¹,⁴*, Sylvan Elhay ², Jochen Deuerlein ³,⁴ and Angus R. Simpson ⁴

¹ INRAE, AQUA Division, UR ETTIS, 50 Av. de Verdun, F-33612 Cestas, France ; olivier.piller@inrae.fr 5

² Sc. of Comp. and Math. Sc., Univ. of Adelaide, Adelaide, SA 5005, Australia; sylvan.elhay@adelaide.edu.au 6

³ 3S Consult GmbH, Garbsen, Germany; deuerlein@3sconsult.de 7

⁴ Sc. of Arch. and Civ. Eng., Univ. of Adelaide, Adelaide, SA 5005, Australia; angus.simpson@adelaide.edu.au 8

* Correspondence: olivier.piller@inrae.fr; 9

† Presented at the 3rd International Joint Conference on Water Distribution Systems Analysis & Computing and Control for the Water Industry (WDSA/CCWI), Ferrara, Italy, July 1st – 4th. 10
11

**Abstract:** First-order approximations have been used with some success for criticality analysis, sensitivity analysis of physical networks, such as water distribution systems, and uncertainty propagation of model parameters. Certain limitations have been reported regarding the accuracy of results, particularly when non-linearity is dominant. In this paper, we show how to efficiently derive the first and second order sensitivities with respect to the variation of their parameters. This makes it possible to improve the first order estimate when necessary. The method is illustrated on a small example system. 12
13
14
15
16
17
18

**Keywords:** Sensitivities; Schur complement; linear equations; sparse matrix; steady state: demand driven modeling; pressure driven modeling; water distribution systems. 19
20
21

## 1. Introduction

Water distribution systems (WDSs) are complex, aging and need to be protected and made more resilient to natural and man-made disasters. There are considerable preservation, health & safety, and sustainability issues at stake in being able to properly manage and understand the operation of such systems. 23
24
25
26

Modeling tools can be very useful for handling such complex systems, to make sustainable management and crisis response decisions. Nevertheless, this may require solving optimization problems using large hydraulic digital models and may prove impossible due to the curse of dimensionality. In response, some authors have suggested 1) using first-order estimates (e.g. the graph Laplacian matrix) [1,2], or 2) using graph partitioning and reduced-order models [3,4] to make the problem tractable. Depending on the problem under consideration, sub-optimality or some kind of limitation may be reported, particularly when precision is required for decision-making and non-linearity is important. 27
28
29
30
31
32
33
34

In addition, uncertainty in the input parameters requires the digital model to be combined with real-time observations to reduce the output uncertainty. Consequently, three main challenges in real-time modeling are 1) reducing computation time, 2) quantifying uncertainties and 3) coupling numerical models with observations. The sensitivity of steady-state solutions to variations in model parameters provides a way of solving the first two challenges [5-7]. 35
36
37
38
39
40

In this research, we show how to derive the first and second order sensitivities of model outputs to variations in parameters by solving linear systems additional to the global gradient algorithm solution. First, we derive explicit formulae for the first and second order sensitivities to parameters. Next, we describe an efficient and low-cost implementation, which uses the Cholesky decomposition of the Schur matrix to calculate 41
42
43
44
45

**Citation:** To be added by editorial staff during production.

Academic Editor: Firstname Last-name

Published: date

<figure><img src="image_1.png" /></figure>

**Copyright:** © 2024 by the authors.
Submitted for possible open access
publication under the terms and
conditions of the Creative Commons
Attribution (CC BY) license
(https://creativecommons.org/licenses/by/4.0/).
sensitivities. Finally, an illustrative example is used to show the potential application for hydraulic modeling of water distribution networks.

## 2. Methods

The method consists of observing that the following general conservation form system can be used to calculate 1st and 2nd order sensitivities of q (the flow rate) and h (the head) wrt parameters:

$$F_{qxy} - Ah_{xy} = x, \quad (1)$$

$$-A^T q_{xy} - Eh_{xy} = y. \quad (2)$$

Where $F = \nabla_q \xi \in \mathbb{R}^{np,np}$ is the Jacobian of the head loss function $\xi$; $A \in \mathbb{R}^{np,nj}$ is the junction node-arc incidence matrix; $E = \nabla_h c \in \mathbb{R}^{nj,nj}$ is the Jacobian of the pressure outflow relationship (POR) function $c$ (for the demand-driven modeling (DDM) case $E = 0$); $x \in \mathbb{R}^{np,nx}$ and $y \in \mathbb{R}^{nj,ny}$ are appropriate vectors or matrices that are specified in Table 1; and $q_{xy}$ (resp.<fnref n="1" /> $h_{xy}$) are flow-rate related quantities (resp. head-related quantities).

<table><caption>**Table 1.** The right-hand sides for system (1-2) and their application.</caption><thead><tr><th>Application</th><th>x</th><th>y</th><th>q<sub>xy</sub></th><th>h<sub>xy</sub></th></tr></thead><tbody><tr><td>Hydraulic state from the linearized system</td><td>e<sup>2</sup></td><td>0</td><td>q<sup>lin</sup></td><td>h<sup>lin</sup></td></tr><tr><td>1<sup>st</sup> order sensitivities wrt demand</td><td>0<sub>np,nj</sub></td><td>∇<sub>d</sub>c</td><td>∇<sub>d</sub>q</td><td>∇<sub>d</sub>h</td></tr><tr><td>1<sup>st</sup> order sensitivities wrt θ, defined in3</td><td>-∇<sub>θ</sub>ξ</td><td>0<sub>nj,ny</sub></td><td>∇<sub>θ</sub>q</td><td>∇<sub>θ</sub>h</td></tr><tr><td>2<sup>nd</sup> order sensitivities wrt demand</td><td>(-∂<sup>2</sup>ξ<sub>j</sub> /∂q<sub>j</sub><sup>2</sup> ∂q<sub>j</sub> /∂d<sub>m</sub> ∂d<sub>n</sub>)</td><td>(∂<sup>2</sup>c<sub>i</sub> /∂h<sub>i</sub><sup>2</sup> ∂h<sub>i</sub> /∂d<sub>m</sub> + E<sub>it</sub> (∂h<sub>i</sub> /∂d<sub>m</sub> + δ<sub>im</sub> ∂h<sub>i</sub> /∂d<sub>n</sub>))</td><td>(∂<sup>2</sup>q<sub>j</sub> /∂d<sub>m</sub> ∂d<sub>n</sub>)</td><td>(∂<sup>2</sup>h<sub>i</sub> /∂d<sub>m</sub> ∂d<sub>n</sub>)</td></tr><tr><td>2<sup>nd</sup> order sensitivities wrt theta</td><td>(-∂<sup>2</sup>ξ<sub>j</sub> /∂q<sub>j</sub><sup>2</sup> ∂q<sub>j</sub> /∂θ<sub>m</sub> ∂θ<sub>n</sub>) - (∂F<sub>ij</sub> /∂r<sub>j</sub>) (∂r<sub>j</sub> /∂θ<sub>n</sub> ∂θ<sub>m</sub> + ∂r<sub>j</sub> /∂θ<sub>m</sub> ∂θ<sub>n</sub>)</td><td>(∂<sup>2</sup>c<sub>i</sub> /∂h<sub>i</sub><sup>2</sup> ∂h<sub>i</sub> /∂θ<sub>m</sub> ∂θ<sub>n</sub>)</td><td>(∂<sup>2</sup>q<sub>j</sub> /∂θ<sub>m</sub> ∂θ<sub>n</sub>)</td><td>(∂<sup>2</sup>h<sub>i</sub> /∂θ<sub>m</sub> ∂θ<sub>n</sub>)</td></tr></tbody></table>

Indeed, system (1-2) is generic as it can be seen in Table 1. if we choose $x = e$ and $y = 0_{nj}$, system (1-2) is the linearized system of pressure-driven modeling (PDM) equations, and $q_{xy}$ and $h_{xy}$ are the estimates of $q$ and $h$ when the head loss and POR models are linear. Likewise, if $x = 0_{np,nj}$ and $y = \nabla_d c$, then $q_{xy} = \nabla_d q$ and $h_{xy} = \nabla_d h$ (see [7] for the derivation). Also, for differentiation wrt $\theta_j = r_j$ or $D_j$ or $\varepsilon_j/D_j$ ($r_j$ the resistance factor, $D_j$ the pipe diameter, $\varepsilon_j/D_j$ the relative roughness of pipe j) the 1st order sensitivities wrt $\theta$ are solutions of (1-2). Just choose $x = -\nabla_\theta \xi$ and $y = 0_{nj,ny}$ (also derived in [7]).

We now consider double scalar differentiation with $d_m$ then $d_n$ (resp. $\theta_m$ then $\theta_n$); the system (1-2) can then be used to calculate the 2nd order sensitivities with appropriate choice of $x$ and $y$ as shown in Table 1. The meaning behind this property is that the flows and heads with their derivatives are sharing similar spatial structures or patterns.

Multiplying Eq. (1) by $A^T F^{-1}$ and adding it to Eq. (2) gives:

$$h_{xy} = -(A^T F^{-1} A + E)^{-1} (A^T F^{-1} x + y), \quad (3)$$

<footnote n="1">¹ resp.: respectively</footnote>

<footnote n="2">² The vector $e$ represents the energy available from source and resource nodes. It is defined as $e = A_0 h_0$.</footnote>

<footnote n="3">³ Where $\theta_j$ is a characteristic parameter of pipe j, such as resistance factor, diameter and relative roughness.</footnote>
It follows from (1) that: 1
$$q_{xy} = F^{-1}(Ah_{xy} + x). \quad (4)$$

Eqs. (3-4) provide a solution template for a linearized estimate of q and h, and the corresponding first order and second order sensitivities. 2
3

Eq. (3) requires the solution of a symmetric matrix equation with the form: 4
$$(A^T F^{-1} A + E) z = w. \quad (5)$$

and this is true for the calculation of all the sensitivities discussed here. If the sensitivities of the solutions to more than one parameter are required, then a significant computational economy can be made. Suppose a first solution is computed using the Cholesky factorization $LL^T = (A^T F^{-1} A + E)$. Sensitivity calculations for any other parameters can be solved with about $2n^2$ floating-point operations each rather than the full Cholesky cost of $O(n^3/6)$ if the same L factor is used with forward- and backward-substitutions. In addition, further savings can be made by exploiting the sparsity of the Cholesky factor. 5
6
7
8
9
10
11

The solution of Eqs (3-4) in this paper was coded with Matlab 2023b. The 1st and 2nd order sensitivities can be calculated for a specific component vector or selected values of interest. This is what we propose in the results. Meanwhile, it is possible to organize the calculation if we are interested in getting an overall view. For example, for the 2nd order sensitivities of q wrt all the demand, there are $nj$ symmetrical matrices $V_{dk}^q$, each of dimension $np \times nj$. Each matrix gives the 2nd order sensitivities of all the flows wrt to all the demands and one single demand. Thus, for example, the matrix for the sensitivities of flows $q_1, q_2, \dots, q_{nj}$ to $d_k$ and all of $d_1, d_2, \dots, d_{nj}$ has the following structure: 12
13
14
15
16
17
18
19

$$V_{d_k}^q = \begin{pmatrix} \frac{\partial^2 q_1}{\partial d_k \partial d_1} & \cdots & \frac{\partial^2 q_1}{\partial d_k \partial d_{k-1}} & \frac{\partial^2 q_1}{\partial d_k^2} & \frac{\partial^2 q_1}{\partial d_k \partial d_{k+1}} & \cdots & \frac{\partial^2 q_1}{\partial d_k \partial d_{nj}} \\ \frac{\partial^2 q_2}{\partial d_k \partial d_1} & \cdots & \frac{\partial^2 q_2}{\partial d_k \partial d_{k-1}} & \frac{\partial^2 q_2}{\partial d_k^2} & \frac{\partial^2 q_2}{\partial d_k \partial d_{k+1}} & \cdots & \frac{\partial^2 q_2}{\partial d_k \partial d_{nj}} \\ \vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots \\ \frac{\partial^2 q_{np}}{\partial d_k \partial d_1} & \cdots & \frac{\partial^2 q_{np}}{\partial d_k \partial d_{k-1}} & \frac{\partial^2 q_{np}}{\partial d_k^2} & \frac{\partial^2 q_{np}}{\partial d_k \partial d_{k+1}} & \cdots & \frac{\partial^2 q_{np}}{\partial d_k \partial d_{nj}} \end{pmatrix}$$

## 3. Results 20

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Example system with two tanks. Pipe lengths are 1000 m; diameters 250 mm; pipe roughness 0.25 mm; source heads are both 100 m; junction node elevations all zero; initial demand at node 1 (resp. node 2) is 60 L/s (resp. 50 L/s). 21
22
23
24</figcaption></figure>

The network used to illustrate the application is shown in **Figure 1**. The Darcy-Weisbach headloss model and the 1-side regularized Wagner model POR with regularization parameter = 1/10 of [7] were used. The demand at node 1 was increased by 5, 10, 20 and 40 L/s. The 2nd order Taylor polynomial approximations to **q** and **h** around the point $d_1$ are given by $X(d_1 + \delta) = X(d_1) + \frac{\partial X}{\partial d_1} \delta + \frac{1}{2!} \frac{\partial^2 X}{\partial d_1^2} \delta^2$. The results for the heads at junction nodes and the flow rates are reported at **Table 2**. We can see the 2nd order estimates are not significantly different from the exact values in column 3. 25
26
27
28
29
30
31
<table><caption>Table 2. 1st and 2nd order estimates for the example network with a demand perturbation of 40 L/s.</caption><thead><tr><th>Name</th><th>h(d<sub>1</sub>)</th><th>h(d<sub>1</sub>+40)</th><th>∂X/∂d<sub>1</sub></th><th>∂<sup>2</sup>X/∂d<sub>1</sub><sup>2</sup></th><th>1<sup>st</sup> order Est.</th><th>Error in m</th><th>2<sup>nd</sup> order Est.</th><th>Error in m</th></tr></thead><tbody><tr><td>Node</td><td>94.64</td><td>89.72</td><td>-0.099382</td><td>-0.001223</td><td>90.67</td><td>-0.95</td><td>89.69</td><td>0.03</td></tr><tr><td>Node</td><td>94.70</td><td>90.71</td><td>-0.089746</td><td>-0.000450</td><td>91.11</td><td>-0.40</td><td>90.75</td><td>-0.04</td></tr><tr><td></td><td>q(d<sub>1</sub>)</td><td>q(d<sub>1</sub>+40)</td><td>∂X/∂d<sub>1</sub></td><td>∂<sup>2</sup>X/∂d<sub>1</sub><sup>2</sup></td><td>1<sup>st</sup> order Est.</td><td>Error in L/s</td><td>2<sup>nd</sup> order Est.</td><td>Error in L/s</td></tr><tr><td>Pipe 1</td><td>55.14</td><td>76.94</td><td>0.524218</td><td>0.001617</td><td>76.11</td><td>0.83</td><td>77.41</td><td>-0.47</td></tr><tr><td>Pipe 2</td><td>-4.86</td><td>-23.06</td><td>-0.475782</td><td>0.001617</td><td>-23.89</td><td>0.83</td><td>-22.59</td><td>-0.47</td></tr><tr><td>Pipe 3</td><td>-54.86</td><td>-73.06</td><td>-0.475782</td><td>0.001617</td><td>-73.89</td><td>0.83</td><td>-72.59</td><td>-0.47</td></tr></tbody></table>

## 4. Discussion and Conclusions

In this paper, the same generic conservative-form system is used to derive linearized estimates of flow and head, and the first order and, for the first time, for the second-order sensitivities. The right-hand-sides of the governing equations change appropriately. Explicit formulae are given and the fact that the same Cholesky factor and sparse solution matrix are shared explains why significant savings can be made in the calculation. It is possible to extend the method to higher order sensitivities.

The development presented in this paper is useful for assessing the probability distributions for link flow rates and nodal piezometric heads. Additionally, it permits Taylor approximation for $q$ and $h$ around known working points. This opens the way to solve difficult problems using a quadratic approximation or to speed up extended period simulations by improving the initial guesses.

**Author Contributions:** Conceptualization and methodology, all authors; software and investigation, S. E., O.P.; validation, S. E.; formal analysis, O.P., S.E., J.D.; original draft preparation, O.P.; writing—review and editing, all; funding acquisition, O.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the French National Research Agency (ANR), grant number ANR-22-CE39-0010, project CoRREau (2024-2027).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The raw data supporting the conclusions of this article will be made available by the authors on request.

**Conflicts of Interest:** The authors declare no conflicts of interest.

## References

1. Abraham, E.; Blokker, E.J.M.; Stoianov, I. Network Analysis, Control Valve Placement and Optimal Control of Flow Velocity for Self-Cleaning Water Distribution Systems. *Procedia Engineering* **2017**, *186*, 576-583, DOI.

2. Steffelbauer, D.B.; Deuerlein, J.; Gilbert, D.; Abraham, E.; Piller, O. Pressure-Leak Duality for Leak Detection and Localization in Water Distribution Systems. *Journal of Water Resources Planning and Management* **2022**, *148*, 04021106, DOI.

3. Elhay, S.; Deuerlein, J.; Piller, O.; Simpson, A.R. Graph Partitioning in the Analysis of Pressure Dependent Water Distribution Systems. *Journal of Water Resources Planning and Management* **2018**, *144*, 04018011-04018011-04018011-04018013, DOI.

4. Pecci, F.; Abraham, E.; Stoianov, I. Model Reduction and Outer Approximation for Optimizing the Placement of Control Valves in Complex Water Networks. *Journal of Water Resources Planning and Management* **2019**, *145*, 04019014, DOI.

5. Fu, G.; Kapelan, Z.; Reed, P. Reducing the Complexity of Multiobjective Water Distribution System Optimization through Global Sensitivity Analysis. *Journal of Water Resources Planning and Management* **2012**, *138*, 196-207, DOI.

6. Hutton, C.; Kapelan, Z.; Vamvakeridou-Lyroudia, L.; Savić, D. Dealing with Uncertainty in WDS Models: A Framework for Real-Time Modeling and Data Assimilation. *Journal of Water Resources Planning and Management* **2014**, *140*, 169-183, DOI.

7. Piller, O.; Elhay, S.; Deuerlein, J.; Simpson, A. Local Sensitivity of Pressure-Driven Modeling and Demand-Driven Modeling Steady-State Solutions to Variations in Parameters. *Journal of Water Resources Planning and Management* **2017**, *143*, 12, DOI.

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.