Dynamical stability of water distribution networks

TL;DR

This paper presents a linear stability analysis of water distribution networks, introducing a resilience index based on eigenvalues that correlates with the system's recovery rate after shocks.

Contribution

It develops a novel analytical framework for assessing the local stability and resilience of nonlinear water distribution network models.

Findings

Steady state is always locally stable.

The resilience index correlates with recovery rate.

Method applicable to arbitrary pipe networks.

Abstract

Water distribution networks are hydraulic infrastructures that aim to meet water demands at their various nodes. Water flows through pipes in the network create nonlinear dynamics on networks. A desirable feature of water distribution networks is high resistance to failures and other shocks given to the system. Such threats would at least transiently change the flow rate in various pipes, potentially undermining the functionality of the whole water distribution system. Here we carry out a linear stability analysis for a nonlinear dynamical system representing the flow rate through pipes that are interconnected through an arbitrary pipe network with reservoirs and consumer nodes. We show that the steady state is always locally stable and develop a method to calculate the eigenvalue that corresponds to the mode that decays the most slowly towards the equilibrium, which we use as an index for resilience of the system. We show that the proposed index is positively correlated with the recovery rate of the pipe network, which was derived from a realistic and industrially popular simulator. The present analytical framework is expected to be useful for deploying tools from nonlinear dynamics and network analysis to designing, resilience managements and scenario testings of water distribution networks.

Equations

$$I_m\frac{\mathrm{d}{Q}_m}{\mathrm{d}t}=h_i-h_j-R_m{Q}_m|Q_m|-u(Q_m).$$

$$R_m=\frac{8f_m{\ell }_m}{g{\pi }^2{d}_m^5},$$

$$f_m={\left(\frac{64}{{\mathcal{R}}_m}\right)}^a{\left(0.75\ln \frac{{\mathcal{R}}_m}{5.37}\right)}^{2(a-1)b}{\left(0.88\ln \frac{6.82d_m}{ϵ}\right)}^{2(a-1)(1-b)}.$$

$${\mathcal{R}}_m=\frac{\left|Q_m\right|d_m}{{\mathcal{A}}_{m}v},$$

$$a=$$

$$b=$$

$$\sum _{j=1;j\ne i}^N{A}_{ij}{Q}_{ij}=-{\tilde{Q}}_i^{\mathrm{ext}}(i=1,2,\dots ,N-N_0),$$

$$\sum _{j=1;j\ne i}^N{A}_{ij}\frac{\mathrm{d}{Q}_{ij}}{\mathrm{d}t}=0(i=1,2,\dots ,N-N_0).$$

$$\sum _{j=1;j\ne i}^N{A}_{ij}\frac{h_i-h_j-R_{ij}{Q}_{ij}|Q_{ij}|-u(Q_{ij})}{I_{ij}}=0(i=1,2,\dots ,N-N_0),$$

$$q_m(Q_m)\equiv R_m{Q}_m|Q_m|+u(Q_m).$$

$$C=\left(\begin{array}{c}{C}^{(1)}\\ C^{(2)}\end{array}\right),$$

$$L_{ij}=\{\begin{array}{ll}{\sum }_{j^{\prime }=1;j^{\prime }\ne i}^N{A}_{i{j}^{\prime }}& (i=j),\\ -A_{ij}& (i\ne j).\end{array}$$

$$L=C{C}^{\top },$$

$$L=\left(\begin{array}{cc}{L}^{(11)}& L^{(12)}\\ L^{(21)}& L^{(22)}\end{array}\right)=\left(\begin{array}{cc}{C}^{(1)}{C}^{(1)\top }& C^{(1)}{C}^{(2)\top }\\ C^{(2)}{C}^{(1)\top }& C^{(2)}{C}^{(2)\top }\end{array}\right),$$

$$L_{\mathrm{w}}=\left(\begin{array}{cc}{L}_{\mathrm{w}}^{(11)}& L_{\mathrm{w}}^{(12)}\\ L_{\mathrm{w}}^{(21)}& L_{\mathrm{w}}^{(22)}\end{array}\right)=CD{C}^{\top }.$$

$$\left(\begin{array}{cc}{L}_{\mathrm{w}}^{(11)}& L_{\mathrm{w}}^{(12)}\end{array}\right)\mathbf{h}-C^{(1)}D\mathbf{q}=0.$$

$$\mathbf{h}=\left(\begin{array}{c}{\mathbf{h}}^{(1)}\\ {\mathbf{h}}^{(2)}\end{array}\right),$$

$$C^{(1)}D{C}^{(1)\top }{\mathbf{h}}^{(1)}=-C^{(1)}D{C}^{(2)\top }{\mathbf{h}}^{(2)}+C^{(1)}D\mathbf{q}.$$

$${\mathbf{h}}^{(1)}={\left(C^{(1)}D{C}^{(1)\top }\right)}^{-1}\left(-C^{(1)}D{C}^{(2)\top }{\mathbf{h}}^{(2)}+C^{(1)}D\mathbf{q}\right).$$

$$\frac{\mathrm{d}\mathbf{Q}}{\mathrm{d}t}=D{C}^{\top }\left(\begin{array}{c}{\left(C^{(1)}D{C}^{(1)\top }\right)}^{-1}\left(-C^{(1)}D{C}^{(2)\top }{\mathbf{h}}^{(2)}+C^{(1)}D\mathbf{q}(\mathbf{Q})\right)\\ {\mathbf{h}}^{(2)}\end{array}\right)-D\mathbf{q}(\mathbf{Q}),$$

$$F_i^{\mathrm{node}}(\mathbf{Q})\equiv \sum _{j=1;j\ne i}^N{A}_{ij}{Q}_{ij}+{\tilde{Q}}_i^{\mathrm{ext}}=0(i=1,2,\dots ,N-N_0).$$

$$F_m^{\mathrm{pipe}}({\mathbf{h}}^{(1)},\mathbf{Q})\equiv h_{i_m}-h_{j_m}-R_m{Q}_m|Q_m|-u(Q_m)=0,$$

$$J=\left(\begin{array}{cc}O& C^{(1)}\\ {\left(C^{(1)}\right)}^{\top }& -\text{diag}(q_m^{\prime }(Q_m))\end{array}\right),$$

$$q_m^{\prime }(Q_m)=2R_m\text{sgn}(Q_m)Q_m+u^{\prime }(Q_m).$$

$$\left(\begin{array}{c}{\mathbf{h}}^{(1)(n+1)}\\ {\mathbf{Q}}^{(n+1)}\end{array}\right)=\left(\begin{array}{c}{\mathbf{h}}^{(1)(n)}\\ {\mathbf{Q}}^{(n)}\end{array}\right)-J^{-1}\left(\begin{array}{c}{\mathbf{F}}_1^{\mathrm{node}}({\mathbf{Q}}^{(n)})\\ ⋮\\ {\mathbf{F}}_{N-N_0}^{\mathrm{node}}({\mathbf{Q}}^{(n)})\\ {\mathbf{F}}_1^{\mathrm{pipe}}({\mathbf{h}}^{(1)(n)},{\mathbf{Q}}^{(n)})\\ ⋮\\ {\mathbf{F}}_M^{\mathrm{pipe}}({\mathbf{h}}^{(1)(n)},{\mathbf{Q}}^{(n)})\end{array}\right).$$

$$J^{\mathrm{dyn}}=\left[D{C}^{(1)\top }{\left(C^{(1)}D{C}^{(1)\top }\right)}^{-1}{C}^{(1)}D-D\right]\text{diag}\left(q_m^{\prime }(Q_m^{\ast })\right).$$

$$B(C^{(1)}{)}^{\top }=0.$$

$${\lambda }_m(B)<{\lambda }_m(B^{\prime }).$$

$${\lambda }_k(J^{\mathrm{dyn}})={\lambda }_k\left({\overline{D}}^{1/2}B{\overline{D}}^{1/2}\right)=$$

$$=$$

$$\le$$

$$=$$

Full text

network, water distribution systems, dynamics, transient flow

I Introduction

Water distribution systems aim to securely provide drinking water to consumers and for firefighting and hence are key infrastructures in our society. They are subjected to drastic changes in demands on the daily and seasonal timescales Butler2006book and to various threats such as contamination Rasekh2014EnvirModelSoft and power outage Gullick2005JWaterSupplyResTechAqua . For safe and secure water supply, it is important that a water distribution system is well designed and operated against potential failures.

In risk management in water engineering, the risk is estimated to be equal to the product of the probability of failures and the consequence of the failure Jun2007JEnvironEng ; Giustolisi2010UrbanWaterJ . In contrast to risk management, resilience management is a relatively new framework to aid us to prepare for potential failures in water distribution systems. Resilience is a broad concept and its definition varies according to studies in engineering Hosseini2016ReliabilityEngSystSafety and in other fields XueWangYang2018NatHazards . The concept of resilience is equally wide for water distribution systems; it is construed as the general capacity of a system to resist, absorb, withstand and recover from stresses Hosseini2016ReliabilityEngSystSafety ; Butler2017GlobalChallenges . It is a useful complement, as many failures that happen in real-life are unforeseeable and are not targeted in the traditional risk management due to the low estimated probability.

Water distribution systems are composed of pipes combined with other functional structures such as reservoirs, pumps and valves. Network disfunctions such as pipe failure are an obvious threat to the functionality of the water distribution system. Network science has accumulated knowledge on how network structure shapes the network’s capacity to withhold random failures and intentional attacks Newman2010book ; Barabasi2016book . Therefore, there has been a surge of interest in water engineering community to deploy network analysis to assess the degree of resilience of water distribution networks. However, to the best of our knowledge, most of the existing proposals of resilience measures for water distribution networks are based on the network structure Yazdani2011Chaos ; Yazdani2011EnvironModelSoftware ; Yazdani2012WaterResourcesRes ; ShuangZhangYuan2014PlosOne ; Herrera2016WaterResourcesMan ; Pundit2016IntJCritInfrastruct ; Hwang2017JWaterResourPlannManag ; Meng2018WaterRes , or steady-state or quasi-steady-state flow rates Todini2000UrbanWater ; Zhuang2013JHydraulEng ; Yannopoulos2013WaterResourManag ; Zhao2015MathProbEng ; Diao2016WaterRes ; Liu2017JWaterResourPlannManag ; as such, the dynamic, transient water flows in the pipe network cannot be represented, which is in fact the key to the resilience performance of a water distribution system. Furthermore, one needs to apply stresses to the system to measure resilience Meng2018WaterRes , which cannot be exhaustive of what failure might occur in a system. Therefore, it is necessary to examine transient dynamics to reveal much of resilience features of the system, in addition to static network structure or flow. In fact, transient dynamics of water flow in pipe networks are not a new issue in water engineering. Examination of transient flows is also a practical requirement because live water distribution networks are incessantly undergoing changes due to, for example, routine operational adjustments, breakdowns of their elements and human error. In turn, transient water flows are known to cause various problems in water distribution systems such as pump failure, collapse of vapor cavities and compromised water quality. Motivated by these needs, various methods to simulate and understand transient flows have been developed Islam1998JHydraulEng ; Boulos2006book ; Todini2011JHydroinfo ; Nault2016JHydraulEng ; JungKarney2017UrbanWaterJ . While useful, these developments are mathematical modelling of the transient water dynamics and their numerical simulation methods. These methods do not themselves provide a measure of resilience.

In the present study, we propose an index of resilience via linear stability analysis of a standard set of dynamical equations modelling the steady-state flow in the given water distribution network. The equations are equivalent to a type of a nonlinear electric circuit system. The proposed resilience index is defined in terms of the dominant eigenvalue of the Jacobian matrix around the steady state. Therefore, the index corresponds to the speed at which the transient dynamics decay towards the steady state in response to a small perturbation in the flow rate. The index is calculated from the structure of the network, the diameter and length of the individual pipes, and other constraints such as the energy level at the reservoir nodes and water demand at individual consumer nodes. We test the proposed index against different resilience measures that we previously proposed based on realistic numerical simulations of water distribution systems. MATLAB code for calculating the so-called local stability index is available at Github (https://github.com/naokimas/LSI_water_distribution_net).

II Model

Consider a connected and undirected network of pipes. By definition, a pipe connects two nodes, where nodes are junctions, reservoirs and households. A dead-end node that is only connected to a single pipe (hence, the node’s degree is equal to one) corresponds to a consumer node such as a household. It should be noted that nodes whose degree is larger than one are junction nodes but can also demand water if they are connected to households. Let us denote by $N$ and $M$ the number of nodes and edges, respectively. Each edge has its own diameter and length, assuming a circular shape of the pipe. These properties affect how much water flows through the pipe under a given condition, thus effectively specifying the weight (i.e., capacity) of the edge. However, we consider the network as unweighted network and explicitly model the effect of the pipe diameter and length.

We consider a rigid water column model on pipe networks under the assumption that flows in a pipe vary slowly in time Wylie1983JHydraulEng ; Islam1998JHydraulEng ; Ghidaoui2005ApplMechRev ; Chaudhry2014book . The dynamical equation representing transient flow through the $m$th pipe ($1\leq m\leq M$), which connects the $i$th and $j$th nodes, is given by

[TABLE]

Nodes $i$ and $j$ depend on edge $m$. However, unless we state otherwise, we use $i$ and $j$ here and in the following text to avoid notational abuse. In Eq. (1), $Q_{m}$ is the time-dependent flow rate through the $m$th pipe, in the direction from the $i$th node to the $j$th node; $t$ represents the time; $h_{i}$ represents the total energy at the $i$th node, which depends on the time; $R_{m}$ is a constant called the pipe coefficient containing the information on the diameter, length and roughness of pipe $m$; $R_{m}Q_{m}|Q_{m}|$ accounts for the nonlinear head loss across the $m$th pipe; $u(Q_{m})$ is the input term which comes from a valve or other structure located on the $m$th pipe; $I_{m}$ is the pipe inertia constant. The dynamical system given by Eq. (1) is equivalent to that of a nonlinear resistor-inductor electrical circuit, where $Q_{m}$ is the current, $h_{i}$ is the voltage, $I_{m}$ is the inductance, $R_{m}$ is a nonlinear resistance and $u(Q_{m})$ is other nonlinear elements placed on the $m$th edge.

The pipe coefficient is given by

[TABLE]

where $f_{m}$ is the dimensionless Darcy-Weisbach friction factor for the $m$th pipe and depends on the flow regime (see below for the formula) and other factors, $\ell_{m}$ is the length of the $m$th pipe, $g$ is the gravitational acceleration and equal to $9.81$ m/s2 and $d_{m}$ is the diameter of the $m$th pipe Boulos2006book . More generally, if the nonlinear head loss term is given by $R_{m}Q_{m}|Q_{m}|^{\gamma-1}$, one obtains $R_{m}=f_{m}\ell_{m}/2gd_{m}\mathcal{A}_{m}^{\gamma}$, where $\mathcal{A}_{m}$ is the cross-sectional area of the $m$th pipe, i.e., $\mathcal{A}_{m}=\pi d_{m}^{2}/4$ Todini1988chapter ; Boulos2006book . The pipe inertia constant is given by $I_{m}=\ell_{m}/(g\mathcal{A}_{m})$.

There are various formula to approximate the friction factor, $f_{m}$, as a function of the pipe diameter, $d_{m}$, the roughness coefficient of the pipe (which is determined by the pipe’s material, age and other factors), Reynolds number and so forth. Because the Reynolds number is a function of the flow rate, $Q_{m}$, the friction factor also depends on $Q_{m}$. Our index developed in section III involves differentiation of the right-hand side of Eq (1) with respect to $Q_{m}$, which requires the derivative of $f_{m}$ with respect to $Q_{m}$. In fact, many formula for $f_{m}$ can be only applicable to particular flow regimes; the flow regime is specified by the value of the Reynolds number. Because $f_{m}$ may depend on pipes, different pipes may be in different flow regimes and some pipes may be situated near the boundary of two flow regimes. Therefore, we select the formula proposed in Ref. Bellos2018JHydraulicEng that is applicable in all flow regimes including laminar and turbulent flows and gives $f_{m}$ as a continuous function of the Reynolds number (and hence continuous in $Q_{m}$). The formula is given by

[TABLE]

In Eq. (3), the Reynolds number for the $m$th pipe is equal to

[TABLE]

where $v=1.007\times 10^{-6}$ m${}^{2}/$s is the kinematic viscosity of water, assuming the temperature of 20 degree celsius. Furthermore, $\epsilon=2.591\cdot 10^{-4}$ m is the roughness coefficient of cast iron Rossman2000EPANET2 (also see Ref. Boulos2006book for a similar value), and

[TABLE]

III Local stability index

In this section, we carry out the linear stability analysis of the steady state of the nonlinear dynamical system given by Eq. (1). The system turns out to be always linearly stable, and we propose a local stability index as the magnitude of the eigenvalue corresponding to the most slowly decaying mode. This task is non-trivial because some nodes act as reservoir of water, other nodes act as consumer of water and the friction factor, $f_{m}$, depends on the dynamical variable, $Q_{m}$.

III.1 Removing redundancy from the dynamical equations

Assume that the total energy at $N_{0}$ nodes is fixed. The $N_{0}$ nodes are typically reservoir whose altitude is higher than typical consumer and junction nodes. Without loss of generality, we assume that $h_{N-N_{0}+1}$, $h_{N-N_{0}+2}$, $\ldots$, $h_{N}$ are fixed. Then, the set of Eq. (1) has $N-N_{0}+M$ unknowns, i.e., $h_{i}$ ($1\leq i\leq N-N_{0}$) and $Q_{m}$ ($1\leq m\leq M$), whereas Eq. (1) only provides $M$ differential equations. The remaining constraints come from the Kirchhoff’s current law (i.e. conservation of water mass) at each of the $N-N_{0}$ nodes whose total energy is not fixed. In other words, we obtain

[TABLE]

where $Q_{ij}$ is the flow rate through pipe $(i,j)$ from the $i$th to $j$th nodes; $\tilde{Q}^{\rm ext}_{i}$ is the external demand (i.e., withdrawal) of water at the $i$th node; $A_{ij}(=A_{ji})$ is the element of the adjacency matrix. In other words, $A_{ij}=1$ if there is a pipe between the $i$th and $j$th nodes. Otherwise, $A_{ij}=0$.

Equations (1) and (7) imply that there are $N-N_{0}+M$ unknowns and the same number of equations for solving the steady state and its linear stability. Now we erase $h_{i}$ ($1\leq i\leq N-N_{0}$) from Eq. (1) using Eq. (7) to derive a self-contained set of $M$-dimensional dynamical system. It should be noted that we do not have to erase $h_{i}$ ($N-N_{0}+1\leq i\leq N$) from Eq. (1) because these $h_{i}$’s are constant. To erase $h_{i}$ ($1\leq i\leq N-N_{0}$), we differentiate Eq. (7) to obtain

[TABLE]

By substituting Eq. (1) in Eq. (8), one obtains

[TABLE]

where $I_{ij}$ is the inertia constant of pipe $(i,j)$. Note that Eq. (2) implies that $R_{ij}=R_{ji}$.

We define the $M$-dimensional diagonal matrix $D$ by the diagonal elements $D_{mm}=1/I_{m}$ ($1\leq m\leq M$). We also denote the head vector by $\bm{h}=(h_{1},\ldots,h_{N})^{\top}$, where $\top$ denotes the transposition. Furthermore, we denote by $\bm{q}=(q_{1},\ldots,q_{M})^{\top}$, of which the $m$th element is given by

[TABLE]

To rewrite Eq. (9), we introduce the incidence matrix, denoted by $C$, which is an $N\times M$ matrix. By definition, corresponding to each $m$th edge $(i,j)$, the entries of the $m$th column of $C$ are given by $C_{i,m}=1$, $C_{j,m}=-1$ and $C_{\ell,m}=0$ ($1\leq\ell\leq N,\ell\neq i,j$). Although the edges are undirected, we assign $+1$ and $-1$ to $i$ and $j$, respectively, by arbitrarily selecting either node forming the $m$th edge to which $+1$ is assigned, for later convenience. We represent $C$ in the block form given by

[TABLE]

where the $(N-N_{0})\times M$ matrix $C^{(1)}$ is the part of the incidence matrix corresponding to the nodes whose total energy is not fixed, and the $N_{0}\times M$ matrix $C^{(2)}$ is the part of the incidence matrix corresponding to the nodes whose total energy is fixed.

We further denote the $N\times N$ combinatorial Laplacian matrix of the network by

[TABLE]

Using

[TABLE]

we obtain the following block format of the Laplacian matrix:

[TABLE]

where $L^{(11)}$ is the $(N-N_{0})\times(N-N_{0})$ matrix corresponding to the nodes whose total energy is not fixed, $L^{(12)}$ is an $(N-N_{0})\times N_{0}$ matrix and so forth.

The combinatorial Laplacian matrix for the network with the same set of edges but with the edge weight given by $D_{mm}=1/I_{m}$ is called the conductance matrix in electrical circuit theory Dorfler2018ProcIeee . We denote the conductance matrix by $L_{\rm w}$. The conductance matrix is given by

[TABLE]

Then, we can rewrite Eq. (9) as

[TABLE]

With the notation

[TABLE]

where $\bm{h^{(1)}}=(h_{1},\ldots,h_{N-N_{0}})^{\top}$ and $\bm{h^{(2)}}=(h_{N-N_{0}+1},\ldots,h_{N})^{\top}$, Eq. (16) is equivalent to

[TABLE]

Matrix $C^{(1)}DC^{(1)\top}$ is so-called the loopy Laplacian matrix Dorfler2018ProcIeee . Because its all eigenvalues are positive if the network is a connected network, which we have assumed, $C^{(1)}DC^{(1)\top}$ has the inverse and one obtains

[TABLE]

This procedure is essentially the same as the network Kron reduction Dorfler2018ProcIeee .

By substituting Eq. (19) in Eq. (1), one obtains

[TABLE]

where $\bm{Q}\equiv(Q_{1},\ldots,Q_{M})^{\top}$. Equation (20) is an $M$-dimensional dynamical system. However, it should be noted that the dynamics are confined on a $(M-N+N_{0})$-dimensional hyperplane specified by Eq. (7), respecting the Kirchhoff’s current law.

III.2 Local stability analysis

In this section, we derive the eigenequation that determines the conventional local linear stability of the steady state of the nonlinear dynamics given by Eq. (1). We use the combination of Eqs. (1) and (7) for solving the steady state and Eq. (20) to perform the local stability analysis around the obtained steady state.

The steady state is given by setting the left-hand side of Eq. (1) to zero. By combining the resulting $M$ constraints with the $N-N_{0}$ constraints provided by Eq. (7), one can obtain the steady state, which consists of the $N-N_{0}+M$ unknowns, i.e., $h_{i}$ ($1\leq i\leq N-N_{0}$) and $Q_{m}$ ($1\leq m\leq M$).

To set up the Newton-Raphson iteration scheme (e.g. Boulos2006book ), we rewrite Eq. (7) as

[TABLE]

We set the left-hand side of Eq. (1) to 0 to obtain

[TABLE]

where $i_{m}$ and $j_{m}$ are the two nodes connected by the $m$th pipe; we use $i_{m}$ and $j_{m}$ instead of $i$ and $j$ to avoid possible notational conflict with Eq. (21). The Jacobian of the right-hand side of Eqs. (21) and (22) with $N+M$ unknowns, $h_{1}$, $\ldots$, $h_{N-N_{0}}$, $Q_{1}$, $\ldots$, $Q_{M}$, is given by

[TABLE]

where $\text{diag}(q_{m}^{\prime}(Q_{m}))$ is the diagonal matrix whose diagonal elements are given by

[TABLE]

The update equation for the Newton-Raphson method is given by

[TABLE]

Denote the obtained steady state by $(h_{1}^{*},\ldots,h_{N-N_{0}}^{*},Q_{1}^{*},\ldots,Q_{M}^{*})$. We now carry out the local stability analysis around this steady state. The Jacobian matrix of the right-hand side of Eq. (20), denoted by $J^{\rm dyn}=(J^{\rm dyn}_{mm^{\prime}})$, is given by

[TABLE]

Proposition III.1**.**

Matrix $B\equiv DC^{(1)\top}\left(C^{(1)}DC^{(1)\top}\right)^{-1}C^{(1)}D-D$ has $(N-N_{0})$-fold zero eigenvalues and $M-N+N_{0}$ negative eigenvalues.

Proof.

A direct substitution verifies

[TABLE]

Therefore, $B$ has $(N-N_{0})$-fold zero eigenvalues, and the corresponding $N-N_{0}$ right eigenvectors are the columns of $(C^{(1)})^{\top}$.

Because $\left(C^{(1)}DC^{(1)\top}\right)^{-1}$ is an $(N-N_{0})\times(N-N_{0})$ matrix, its rank is at most $N-N_{0}$. Therefore, the $M\times M$ matrix $B^{\prime}\equiv C^{(1)\top}\left(C^{(1)}DC^{(1)\top}\right)^{-1}C^{(1)}D$ has at most rank $N-N_{0}$ and therefore has at least $M-N+N_{0}$ zero eigenvalues.

We denote the eigenvalues of a symmetric matrix $X$ by $\lambda_{1}(X)\leq\cdots\leq\lambda_{M}(X)$. Using the fact that $D$ is positive definite because all its diagonal elements are positive (i.e., $=1/I_{m}$), the Weyl’s theorem on eigenvalues applied to the sum of $B^{\prime}$ and $-D$, which is equal to $B$, implies that

[TABLE]

Note that $\lambda_{m_{0}}(B^{\prime})=\lambda_{m_{0}+1}(B^{\prime})=\cdots=\lambda_{m_{0}+M-N+N_{0}-1}(B^{\prime})=0$ for an $m_{0}$ ($1\leq m_{0}\leq N-N_{0}+1$). Therefore, Eq. (28) yields $\lambda_{m_{0}}(B),\ldots,\lambda_{m_{0}+M-N+N_{0}-1}(B)<0$. Because this relationship must be consistent with the fact that $B$ has $(N-N_{0})$-fold zero eigenvalues, we obtain $m_{0}=1$ and $\lambda_{1}(B),\ldots,\lambda_{M-N+N_{0}}(B)<0,\lambda_{M-N+N_{0}+1}(B)=\cdots=\lambda_{M}(B)=0$. ∎

Let us assume for any pipe $m$ that $q_{m}^{\prime}(Q_{m}^{*})>0$, which means that the head loss owing to friction along the pipe increases as the steady-state flow rate, $Q_{m}^{*}$, increases. We verified $q_{m}^{\prime}(Q_{m}^{*})>0$ held true in all the numerical simulations carried out in the next section. Then, $J^{\rm dyn}$ has $M-N+N_{0}$ negative eigenvalues and $N-N_{0}$ zero eigenvalues as matrix $B$ does (Appendix A).

Therefore, the dynamics given by Eq. (20) are locally neutrally stable in the subspace spanned by the column vectors of $\left[\text{diag}\left(q_{m}^{\prime}(Q_{m}^{*})\right)\right]^{-1}(C^{(1)})^{\top}$. However, perturbations along these directions are infeasible because such a perturbation is not on the hyperplane defined by Eq. (7) on which the dynamics are constrained. Note that the normal vectors of the hyperplane are given by the row vectors of $C^{(1)}$, i.e., column vectors of $(C^{(1)})^{\top}$. Therefore, in the local stability analysis, we should ignore these $N-N_{0}$ zero eigenvalues and inspect the other $M-N+N_{0}$ eigenvalues. With this consideration, we conclude that the steady state constrained by the Kirchhoff’s current law (i.e., Eq. (7)) is always locally stable.

To define a local stability index, we calculate the largest eigenvalue of $J^{\rm dyn}$ except the $(N-N_{0})$-fold zero eigenvalue. Note that all eigenvalues of $J^{\rm dyn}$ are real although $J^{\rm dyn}$ is not a symmetric matrix in general (Appendix A). We consider the the largest non-zero eigenvalue of $J^{\rm dyn}$ (which is negative) and refer to its absolute value as the local stability index, denoted by $\rho$. A large $\rho$ value implies that dynamics perturbed by a small amount from the steady state would rapidly return to the steady state. Mathematically, the perturbations are given to the flow rate for some pipes. However, we expect that $\rho$ serves as an index that can characterise dynamical stability of water distribution networks in wider contexts.

IV Numerical results

In this section, we calculate the local stability index for various networks and examine its relevance to stability and resilience in simulated dynamics of water flow.

We use the synthesized networks used in our previous study Meng2018WaterRes . As supplementary materials of Ref. Meng2018WaterRes , the full network structure and properties of system components (e.g., pipes, nodal demands, elevations) were made open to the public in the EPANET2 format, except for the one network. For completeness, we have made the complete data for all the networks and their system components used in the present study available in the MATLAB format at Github, alongside the code for calculating the local stability index.

Of the 85 networks used in our previous study Meng2018WaterRes , five networks include pumps, which would pin the flow rate to a constant value. Because our formalism does not directly cover the case of the pinned flow rate, we exclude these five networks. Of the remaining 80 networks, 10 networks have $N=102$ nodes, $N_{0}=2$ reservoir nodes among the $N$ nodes and $M=110$ pipes; 11 networks have $(N,N_{0},M)=(102,2,130)$; 10 networks have $(N,N_{0},M)=(204,4,223)$; 10 networks have $(N,N_{0},M)=(204,4,263)$; 9 networks have $(N,N_{0},M)=(306,6,335)$; 9 networks have $(N,N_{0},M)=(306,6,395)$; 5 networks have $(N,N_{0},M)=(404,4,443)$; 5 networks have $(N,N_{0},M)=(404,4,523)$; 5 networks have $(N,N_{0},M)=(406,6,525)$; 5 networks have $(N,N_{0},M)=(408,8,447)$; 1 network has $(N,N_{0},M)=(506,6,554)$. We calculated $\rho$ for each of the 80 networks.

The diameter of a majority of pipes was equal to 400 mm. There were also a few larger values of the diameter (e.g., 900 mm). Such wider pipes were incident (i.e., directly connected) to a reservoir. The diameter of the pipes was automatically generated by the software Decorte2014WaterResourcesMan .

We set the total energy of all the $N_{0}$ reservoirs in each network to 65 m Decorte2014WaterResourcesMan . In addition to these $N_{0}$ nodes, some nodes in each network had zero water demand (i.e., $\tilde{Q}_{i}^{\rm ext}=0$) because they are pure junctions connecting different pipes. We set $u(Q_{m})$ for each pipe to zero because the data are from the 80 networks that do not have valves or pumps.

In our previous study, we numerically simulated dynamics of water flow in these networks using software EPANET2 Mugume2015WaterRes ; Diao2016WaterRes and measured six strain indices to characterise resilience of water distribution systems. These indices are as follows (see Meng2018WaterRes for fuller definitions). First, the time to strain is defined as the time between the application of stress and the start of service failure, which is when the level of service at a node drops below a predefined threshold. Second, the failure duration is defined as the time needed for the system to recover to normal performance since the service failure triggered by the application of stress. Third, the failure magnitude is defined as the most severe drop at a node in terms of the system service performance as a result of the administered stress, averaged over all the nodes. Fourth, the failure rate is defined as the failure magnitude divided by the time between the start of the failure and the occurrence of the worst system performance. Fifth, the recovery rate is defined as the failure magnitude divided by the time between the worst system performance and the return to the performance threshold value used in the definition of the time to strain. Sixth, the severity is defined as the threshold minus the system performance (which is positive during the system failure) integrated with respect to time over the period of system failure.

The relationship between each of the six strain indices and the local stability index, $\rho$, is shown in Fig. 1. In the figure, each circle represents one of the 80 networks, $r$ represents the Pearson correlation coefficient and $p$ represents the $p$-value for the Pearson correlation coefficient. We found a positive correlation between $\rho$ and the recovery rate (Fig. 1(e)). This result is intuitive because the recovery rate quantifies how the system returns to the normal performance after being perturbed by external stress, which agrees with the concept of linear stability when the given external stress is infinitesimally small. The obtained correlation is only moderate (i.e., $r=0.420$, $p\approx 1.04\times 10^{-4}$) probably due to various nonlinearity of the system, which our linear stability analysis does not aim to explain. However, we consider that $\rho$ captures a strain index derived from EPANET2, which is an involved numerical simulation package, reasonably well.

Figure 1 also indicates that $\rho$ is strongly correlated with the failure magnitude (Fig. 1(c)) and the failure rate (Fig. 1(d)). We do not have an active interpretation of this result. In fact, it apparently sounds contradictory that a large $\rho$ value, which predicts that the network is resilient, is associated with a large failure magnitude and failure rate. However, the failure magnitude and failure rate quantify the magnitude of failure response induced by external stress, which does not have to do with our local stability analysis even in the limit of infinitesimally small external stress. Therefore, we conclude that the failure magnitude and the failure rate characterise different aspects of the system’s resilience from what the recovery rate and $\rho$ do. This view is consistent with our additional observation that, within our previous numerical results Meng2018WaterRes , there is positive correlation between the failure magnitude and the recovery rate ($r=0.917$, $p<10^{-32}$) and between the failure rate and the recovery rate ($r=0.677$, $p\approx 5.60\times 10^{-12}$).

In our previous study, we measured eight structural properties of the water distribution network and examined the association between each of them and each of the six strain indices Meng2018WaterRes . Here we measured the Pearson correlation between each of these structural properties, whose definitions are in Appendix B, and the recovery rate. We found that none of the eight structural properties of the network was as strongly correlated with the recovery rate as $\rho$ was (link density: $r=0.145$, $p=0.204$; algebraic connectivity: $r=0.114$, $p=0.316$; diameter: $r=-0.143$, $p=0.207$; average path length: $r=-0.204$, $p=0.0717$; central point dominance: $r=-0.0515$, $p=0.652$; heterogeneity: $r=0.335$, $p=0.0026$; spectral gap: $r=-0.0509$, $p=0.656$; clustering coefficient: $r=0.198$, $p=0.0798$). Therefore, the local stability index is better at capturing the recovery rate than these structural properties. This is probably because these structural properties do not directly quantify the speed at which the water distribution network responds to a perturbation to the steady state.

V Discussion

We carried out a local stability analysis of water flow in pipe networks including reservoir nodes and consumer nodes. The steady state was shown to be always linearly stable. We used the eigenvalue corresponding to the slowest relaxation mode as the local stability index. The proposed index was moderately correlated with the recovery rate in response to external stress, which was numerically obtained from an involved simulator commonly used in water engineering research community.

We used a particular formula for the friction factor Bellos2018JHydraulicEng . However, the present framework accepts other formulae for the friction factor. Because the differentiation of the friction factor with respect to the Reynolds number is used in calculating the local stability index and the value of the friction factor depends on the pipe in general, one should use a formula that is continuous and covers a wide range of the Reynolds number and hence different flow schemes. Some possible choices of the formula different from the one used in the present paper are found in the literature Swamee1993JTransportEng ; Cheng2008JHydraulicEng ; Brkic2018ApplSci . It is straightforward to extend our local stability index to the case of different functional forms of the friction factor.

In 1972, May proposed to examine the eigenvalue spectra of interaction networks to characterise dynamical stability of complex systems, in particular ecological systems May1972Nature-stability . This is a landmark example of local stability analysis. While the original analysis and many studies that ensued were for random networks May1972Nature-stability , the understanding has been extended to the case of interaction networks with general network structure Okuyama2008EcolLett ; Allesina2012Nature . This family of method is similar to the one proposed in the present study. Main differences are that we have started from a particular set of hydraulic equations modelling water flow in pipe networks and that we have found that our system is always locally stable and hence used the relaxation mode to quantify the extent of the local stability of the system.

The proposed local stability index was significantly correlated with the recovery rate, which was obtained from numerical simulations in which the applied shock was not necessarily small. However, because the present analysis is a local stability analysis, it does not tell us the behaviour of the system when a shock injected to the steady state is not small. In water engineering applications, the magnitude of perturbation is not necessarily small Gullick2005JWaterSupplyResTechAqua . A recent study in network science has developed a formalism that reduces a class of dynamical systems on interaction networks into a one-dimensional effective dynamical system and assesses the resilience of the dynamical system against perturbation that is not necessarily small GaoBarzelBarabasi2016Nature . Therefore, adapting their resilience formalism to the case of transient and steady-state dynamics of water flow in pipe networks is of a practical relevance. In Eq. (1), there are usually only a fraction of nodes whose energy value is fixed (i.e., reservoirs) and different nodes have different water demands. In this sense, the flow dynamics on pipe networks are heterogeneous in addition to being heterogeneous in the network structure and edge weight. Other structures such as valves would add more complexity and heterogeneity to the system. For such a heterogeneous system, it is unclear whether one can transform Eq. (1) into a form that is compatible with the resilience function formalism GaoBarzelBarabasi2016Nature ; Tu2017PhysRevE . This issue is left as future work.

Data accessibility

The data on the network structure and system components are available on Github.

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

NM and FM conceived of and designed the study. FM provided the data on the network structure and system components. NM carried out the theoretical analysis and performed the data analysis. NM and FM discussed the results and drafted the manuscript. All authors read and approved the manuscript.

Funding

N.M. and F.M. acknowledge the support provided through BRIM (Building Resilience into Risk Management) project, EPSRC (EP/N010329/1).

Appendix A: Eigenvalues of $J^{\rm dyn}$

Here we provide an elementary proof that $J^{\rm dyn}$ has $M-N+N_{0}$ zero eigenvalues and $N-N_{0}$ negative eigenvalues when $q_{m}^{\prime}(Q_{m}^{*})>0$, $1\leq m\leq M$.

Let us introduce a short-hand notation $\overline{D}\equiv\text{diag}\left(q_{m}^{\prime}(Q_{m}^{*})\right)$. Then, one obtains $J^{\rm dyn}=B\overline{D}$. Because $J^{\rm dyn}$ is similar to the symmetric matrix $\overline{D}^{1/2}B\overline{D}^{1/2}$, the min-max theorem for the $k$th smallest eigenvalue implies that

[TABLE]

Similarly, one obtains

[TABLE]

By combining $\lambda_{1}\left(\overline{D}\right)>0$, $\lambda_{M}\left(\overline{D}\right)>0$ and $\lambda_{1}(B),\ldots,\lambda_{M-N+N_{0}}<0,\lambda_{M-N+N_{0}+1}=\cdots=\lambda_{M}=0$, one obtains $\lambda_{1}(J^{\rm dyn}),\ldots,\lambda_{M-N+N_{0}}(J^{\rm dyn})<0,\lambda_{M-N+N_{0}+1}(J^{\rm dyn})=\cdots=\lambda_{M}(J^{\rm dyn})=0$.

Appendix B: Structural properties of networks

The eight structural properties of the network that we measured in our previous study Meng2018WaterRes and used in the present article for comparison purposes are as follows.

The link density is given by $2M/\left[N(N-1)\right]$. The algebraic connectivity is the smallest positive eigenvalue of the Laplacian matrix, $L$, of the undirected and unweighted network. The version of the clustering coefficient used in Meng2018WaterRes is given by $3\times(\text{number of triangles in the network})/(\text{number of connected triples of nodes in the network})$. The average path length is equal to the shortest path length between two nodes, which is averaged over all the possible $N(N-1)/2$ node pairs. The central point dominance is given by $\sum_{i=1}^{N}(b_{\max}-b_{i})/(N-1)$, where $b_{i}$ is the betweenness centrality of the $i$th node and $b_{\max}=\max_{j=1,\ldots,N}b_{j}$. The heterogeneity is that in the degree and given by $\sum_{i=1}^{N}\left(k_{i}-\langle k\rangle\right)^{2}/\langle k\rangle$, where $k_{i}$ is the degree of the $i$th node and $\langle k\rangle=\sum_{i=1}^{N}k_{i}/N$ is the average degree of the node. The version of the spectral gap used in Meng2018WaterRes is given by the difference between the two largest eigenvalues of $L$. The modularity is a quantity representing the quality of community structure detected in the network and approximate modularity maximisation was carried out using Newman’s algorithm Newman2004PhysRevE-fast .