On the Flow Problem in Water Distribution Networks: Uniqueness and Solvers

TL;DR

This paper investigates the water flow problem in distribution networks, establishing solution uniqueness, extending convex relaxation methods, and proposing multiple solvers with guaranteed convergence for different network configurations.

Contribution

It introduces new theoretical results on solution uniqueness and develops a hierarchy of solvers, including convex relaxations and algorithms, for water flow problems in complex water distribution networks.

Findings

Uniqueness of water flow solutions is extended to networks with multiple fixed-pressure nodes.

Convex relaxation methods are effective for networks without cycles and with certain pump configurations.

The proposed solvers reliably converge regardless of initial conditions, validated on benchmark networks.

Abstract

Increasing concerns on the security and quality of water distribution systems (WDS), call for computational tools with performance guarantees. To this end, this work revisits the physical laws governing water flow and provides a hierarchy of solvers of complementary value. Given the water injection or pressure at each WDS node, finding the water flows within pipes and pumps along with the pressures at all WDS nodes constitutes the water flow (WF) problem. The latter entails solving a set of (non)-linear equations. We extend uniqueness claims on the solution to the WF equations in setups with multiple fixed-pressure nodes and detailed pump models. For networks without pumps, the WF solution is already known to be the minimizer of a convex function. The latter approach is extended to networks with pumps but not in cycles, through a stitching algorithm. For networks with non-overlapping cycles, a provably exact convex relaxation of the pressure drop equations yields a mixed-integer quadratically-constrained quadratic program (MI-QCQP) solver. A hybrid scheme combining the MI-QCQP with the stitching algorithm can handle WDS with overlapping cycles, but without pumps on them. Each solver is guaranteed to converge regardless of initialization, as numerically validated on a benchmark WDS.

Figures

Equations

$$d_n=\sum _{k:(n,k)\in \mathcal{P}}{f}_{nk}-\sum _{k:(k,n)\in \mathcal{P}}{f}_{kn}.$$

$$A_{p,k}=\{\begin{array}{ll}+1,& k=m\\ -1,& k=n\\ 0,& \text{otherwise}\end{array}\forall p=(m,n)\in \mathcal{P}.$$

$${\mathbf{A}}^{\top }\mathbf{f}=\mathbf{d}.$$

$$h_m-h_n=c_{mn}\mathrm{sign}(f_{mn})f_{mn}^2$$

$$\operatorname{sign}(x):=\left\{\begin{array}[]{rl}+1,&~{}x>0\\ 0,&~{}x=0\\ -1,&~{}x<0\end{array}\right..$$

$$h_n-h_m=g_{mn}.$$

$$g_{mn}(f_{mn},{\omega }_{mn})={\lambda }_{mn}{f}_{mn}^2+{\mu }_{mn}{\omega }_{mn}{f}_{mn}+{\nu }_{mn}{\omega }_{mn}^2$$

$$g_{mn}(f_{mn};{\omega }_{mn}^0)=h_n-h_m={\lambda }_{mn}{f}_{mn}^2+{\bar{\mu }}_{mn}{f}_{mn}+{\bar{\nu }}_{mn}$$

$$\mathrm{find}$$

$$s.to$$

$${\mathbf{n}}^{\top }\mathbf{A}(\tilde{\mathbf{h}}-\mathbf{h})=0.$$

$${\mathbf{n}}_{+}^{\top }\left(\mathbf{A}\tilde{\mathbf{h}}-\mathbf{Ah}\right)={\mathbf{n}}_{-}^{\top }\left(\mathbf{A}\tilde{\mathbf{h}}-\mathbf{Ah}\right).$$

$${\mathbf{a}}_p^{\top }\tilde{\mathbf{h}}>{\mathbf{a}}_p^{\top }\mathbf{h}\text{if and only if}{\tilde{f}}_p>f_p$$

$${\mathbf{n}}_{+}^{\top }\left(\mathbf{A}\tilde{\mathbf{h}}-\mathbf{Ah}\right)>0.$$

$${\mathbf{n}}_{-}^{\top }\left(\mathbf{A}\tilde{\mathbf{h}}-\mathbf{Ah}\right)<0.$$

$$\underset{\mathbf{f}}{\min }$$

$$s.to$$

$$\underset{\{f_p,w_p,y_p,t_p\}}{\min }$$

$$s.to$$

$$-w_p\le f_p\le w_p,$$

$$w_p^2\le y_p,y_p^2\le w_p{t}_p,$$

$$f_p^{k+1}:=\arg \underset{f_p}{\min }\frac{c_p|f_p{|}^3}{3}-f_p{\mathbf{a}}_p^{\top }{\mathbf{\xi }}^k=\mathrm{sign}({\mathbf{a}}_p^{\top }{\mathbf{\xi }}^k)\sqrt{\frac{|{\mathbf{a}}_p^{\top }{\mathbf{\xi }}^k|}{c_p}}.$$

$${\mathbf{\xi }}^{k+1}:={\mathbf{\xi }}^k+\mu \left(\mathbf{d}-{\mathbf{A}}^{\top }{\mathbf{f}}^{k+1}\right).$$

$$\underset{\mathbf{h}}{\min }\frac{2}{3}\sum _{(m,n)\in \mathcal{P}}\frac{|h_m-h_n{|}^{\frac{3}{2}}}{\sqrt{c_{mn}}}-{\mathbf{d}}^{\top }\mathbf{h}.$$

$$[{\mathbf{g}}_u(\mathbf{h}){]}_n=\sum _{p=(m,n)\in \mathcal{P}}\mathrm{sign}({\mathbf{a}}_p^{\top }\mathbf{h})\sqrt{\frac{|{\mathbf{a}}_p^{\top }\mathbf{h}|}{c_p}}-d_n.$$

$$\underset{\mathbf{f}}{\min }$$

$$s.to$$

$$\underset{\mathbf{h}}{\min }\frac{\rho }{\rho +1}\sum _{(m,n)\in \mathcal{P}}\frac{|h_m-h_n{|}^{\frac{\rho +1}{\rho }}}{c_{mn}^{1/\rho }}-{\mathbf{d}}^{\top }\mathbf{h}.$$

$$-$$

$$h_m-h_n+c_{mn}{f}_{mn}^2\le M{x}_{mn}$$

$$x_{mn}\in \{0,1\}.$$

$$h_n-h_m\le {\lambda }_{mn}{f}_{mn}^2+{\hat{\mu }}_{mn}{f}_{mn}+\hat{\nu }$$

$$\underset{\mathbf{f},\mathbf{h},\mathbf{x}}{\min }$$

$$s.to$$

Full text

I Introduction

Water distribution systems serve as a critical infrastructure across the world. The direct dependence of human lives on the availability of water has motivated research on the security, resiliency, and quality of water supply systems [1], [2]. The high cost of installation for different WDS components renders long-term network planning an important problem [3], [4]. Furthermore, the relatively expensive operation of a WDS is primarily attributed to the electricity cost for running pumps to properly circulate water [5]. Thus, optimal pump scheduling for the daily operation of a WDS is a pertinent research problem [6], [7], [8], [9]. An inevitable component of the aforementioned computational problems is satisfying the physical laws governing water flow. Mathematically, the water flow (WF) equations consist of a set of linear equations ensuring mass conservation, along with a set of non-linear equations arising from energy and momentum conservation [10].

Solving the WF equations constitutes the water flow problem. Specifically, given water demand at all nodes, the standard WF task aims at finding the water flows within all pipes and the pressures at all nodes complying with the WF equations. Despite nonlinear, these equations enjoy a unique solution for networks with a single fixed-pressure node and when the pressure added by pumps is approximated as constant [11]. Modern renditions of the WF problem may incorporate pumps, valves, and pressure-based demands [12], [13]. Either way, handling the non-linear equations stemming from the conservation of energy and momentum remain the core challenge [10].

Existing WF solvers update iteratively a set of WF variables, which could be the pipe flows, loop flows, nodal pressures, or combinations thereof [14]. These solvers can be broadly classified into those relying on successive linear approximations, and those relying on Newton-Raphson-type of updates [10]. For example, the WF solver of [15] constitutes a fixed-point iteration and belongs to the former class, while EPANET (perhaps the most widely used WF solver) to the latter [13]. In fact, most of the schemes within each class have been shown to be equivalent to each other upon a (non)-linear transformation of variables [14], [10]. To cope with the problem of dimensionality, different preprocessing, partitioning, and reformulations have been reported in [16], [17], [18].

The aforesaid solvers exhibit two major shortcomings. First, their convergence and rate of convergence depend critically on initialization. For example, it has been numerically demonstrated that EPANET fails to find a WF solution for some practical water networks [19], [20]. Nonetheless, proper initialization may be challenging when dealing with stochastic planning or risk analysis, where the WF task has to be solved repeatedly and under varying demands [17]. As a second shortcoming, the existing solvers do not naturally extend to optimal water flow (OWF) formulations. Therefore, most OWF efforts resort to linear approximations; non-linear local optimization; or slow zero-order algorithms building on an independent WF solver such as EPANET; see [1]. Such OWF approaches lack scalability and/or optimality guarantees.

The contribution of this work is in two fronts. After a brief modeling of water networks (Section II), this work first generalizes the uniqueness of the WF solution to WDS setups with multiple fixed-pressure nodes and more practical pump models with flow-dependent pressure gains. Second, it puts forth a suite of solvers that can provably recover the WF solution under different network setups; see also Fig. 1:

1. i)

In networks without pumps, it is already known that the WF solution can be recovered as the minimizer of a convex problem [11]. Sections IV-A devises SOCP- and dual decomposition-based solvers, while Section IV-C applies them to other pressure drop laws. 2. ii)

This energy function-based approach is extended to networks where pumps do not lie on cycles through the stitching algorithm of Section IV-B. 3. iii)

In networks with no overlapping cycles, the mixed-integer quadratically-constrained quadratic program (MI-QCQP) of Section V can find the WF solution. 4. iv)

A hybrid procedure combining the stitching algorithm and the MI-QCQP solver handles networks having overlapping cycles, but without pumps on them (Section VI).

These novel solvers not only operate under different network configurations, but constitute a hierarchy: Solver ii) builds on i); and solver iv) builds on ii) and iii). Numerical tests on benchmark networks evaluate the correctness and running times for i) and iii); and validate them against the EPANET solver (Section VII).

Regarding notation, lower- (upper-) case boldface letters denote column vectors (matrices). Calligraphic symbols are reserved for sets. The vectors of all zeros, all ones, and the $n$-th canonical vector are denoted respectively by $\mathbf{0}$, $\mathbf{1}$, and $\mathbf{e}_{n}$, while their dimension will be clear from the context. The symbol ⊤ stands for transposition.

II Water Distribution System Modeling

A WDS can be represented by a directed graph ${\mathcal{G}:=(\mathcal{N},\mathcal{P})}$. Its nodes are indexed by $n\in\mathcal{N}:=\{1,\dots,N\}$ and correspond to water reservoirs, tanks, and points of water demand. Let $d_{n}$ be the rate of water injected into the WDS from node $n$. For reservoirs apparently $d_{n}\geq 0$; for nodes with water consumers $d_{n}\leq 0$; tanks may be filling or emptying; and $d_{n}=0$ for junction nodes.

The edges in set $\mathcal{P}$ with cardinality $P:=|\mathcal{P}|$, are associated with pipes and pumps. The directed edge $p=(m,n)\in\mathcal{P}$ models the water pipe between nodes $m$ and $n$. Its water flow is denoted by $f_{mn}$ or $f_{p}$ depending on the context. If water flows from node $m$ to $n$, then $f_{mn}\geq 0$; otherwise $f_{mn}<0$.

Conservation of water flow dictates that for all $n\in\mathcal{N}$

[TABLE]

The connectivity of the WDS is captured by the edge-node incidence matrix $\mathbf{A}\in\mathbb{R}^{P\times N}$ with entries

[TABLE]

Given $\mathbf{A}$ and upon stacking flows and injections respectively in $\mathbf{f}\in\mathbb{R}^{P}$ and $\mathbf{d}\in\mathbb{R}^{N}$, the conservation of water flow across the WDS can be compactly expressed as

[TABLE]

The operation of WDS is also governed by pressures. Water pressure is surrogated by pressure head, defined as the equivalent height of a water column in meters, which exerts the surrogated pressure at its bottom. The pressure heads at all WDS nodes are measured with respect to a common geographical elevation level. The pressure head (henceforth pressure) at node $n$ is denoted by $h_{n}$. Moreover, the pressure $h_{r}$ at a reservoir node $r\in\mathcal{N}$ typically serves as the pressure of reference.

With water flowing in a pipe, pressure drops along the direction of flow due to friction. The pressure drop across pipe $(m,n)\in\mathcal{P}$ is described by the Darcy-Weisbach law [21]

[TABLE]

where the constant $c_{mn}$ depends on pipe dimensions [8], and the sign function is defined as

[TABLE]

The pressures at all nodes are collected in vector $\mathbf{h}\in\mathbb{R}^{N}$.

To maintain pressures at desirable levels, water utilities use pumps on specific pipes. Let $\mathcal{P}_{a}\subset\mathcal{P}$ be the subset of pipes hosting a pump. The pipes in $\mathcal{P}_{a}$ can be considered lossless; this is without loss of generality since a pump can be modeled by an ideal pump followed by a short pipe. The remaining edges form the subset $\bar{\mathcal{P}}_{a}:=\mathcal{P}\setminus\mathcal{P}_{a}$, and correspond to lossy pipes governed by (3). When pump $p=(m,n)\in\mathcal{P}_{a}$ is running, it adds pressure $g_{mn}\geq 0$ so that

[TABLE]

As detailed below, the pressure added by a pump decreases with the water flow through the pump [22], [23]. Moreover, for variable speed pumps, the pressure added increases with the pump speed. The exact relation is provided by manufacturers in the form of pump operation curves, and are oftentimes approximated with quadratic curve fits [22], [23], [7]. In detail, the pressure added by pump $p=(m,n)$ is modeled as

[TABLE]

where $\omega_{mn}$ is the pump speed; and $\lambda_{mn}<0$, $\mu_{mn}\geq 0$, and $\nu_{mn}\geq 0$ are known pump parameters [7]. When pump $(m,n)$ is running, its flow has to be maintained within the range $0\leq\underline{f}_{mn}\leq f_{mn}\leq\bar{f}_{mn}$ due to engineering limitations [21]. Whereas the speed for fixed-speed pumps is constant; for variable-speed pumps, it is controlled by the operator. Either way, for a WF task the speed of every pump is given; otherwise the problem would become under-determined. Given its speed $\omega_{mn}^{0}$, the pressure added by pump $(m,n)$ is

[TABLE]

where $\bar{\mu}_{mn}:=\mu_{mn}\omega_{mn}^{0}$ and $\bar{\nu}_{mn}:=\nu_{mn}(\omega_{mn}^{0})^{2}$. Because the pump parameters satisfy $2\lambda_{mn}\underline{f}_{mn}+\mu_{mn}\omega_{mn}<0$ over the operating range of pump speeds $\omega_{mn}$, the pressure gain due to any pump is a strictly decreasing function of the water flow in the range $[\underline{f}_{mn},\bar{f}_{mn}]$. This observation is instrumental in establishing the uniqueness of the WF solution in Section III.

When pump $p=(m,n)\in\mathcal{P}_{a}$ is not running, water can flow freely in either directions through a bypass valve [22]. Because bypass valve sections are typically short, one can ignore the pressure drop along them to get $h_{m}=h_{n}$ and $g_{mn}=0$. In this case, the WDS graph can be reduced by removing pipe $p$ and node $n$, and connecting to node $m$ the edges previously incident to $n$. Alternatively, the valve can be modeled as a short lossy pipe, whose pressure drop is governed by (3).

Valves constitute a vital component for water flow control. They can be modeled by an on/off switch; a linear pressure-reducing model; a flow-dependent non-linear model; or a flow control model [7], [13]. Under the typical operational setup, the valve at the reference node regulates pressure, whereas the valves at the remaining reservoirs and tanks regulate flows.

Summarizing this section, a WDS operating point is described by the triplet $(\mathbf{d},\mathbf{f},\mathbf{h})$ satisfying the WF equations of (2), (3), and (5). This work deals with the uniqueness of a WF solution and efficient solvers for finding this solution.

III Uniqueness of Water Flow Solution

Different from the OWF problem where tanks and pumps are scheduled over a time horizon, the WF task aims at solving the WF equations given the water injections $\mathbf{d}$. As such, it constitutes a key component of WDS operation and planning. The WF problem is formally stated next.

Definition 1**.**

Given: i) water injections $\mathbf{d}$; ii) the statuses and speeds $\{\omega_{mn}\}$ for all pumps $(m,n)\in\mathcal{P}_{a}$; and iii) the pressure $h_{r}$ at the reference node $r\in\mathcal{N}$; the WF task aims at finding the flows $\mathbf{f}$ and pressures $\mathbf{h}$ satisfying (2), (3), (5).

The WF task involves $N+P-1$ equations over $N+P-1$ unknowns. The water balance in (2) yields $N-1$ linearly independent equations. In addition, the pressure drops across lossy pipes [cf. (3)], and the pressure gains due to pumps [cf. (5)] provide $P$ non-linear equations.

According to Definition 1, both $\mathbf{f}$ and $\mathbf{h}$ are unknown. Nevertheless, to find a triplet $(\mathbf{d},\mathbf{f},\mathbf{h})$ that satisfies (2), (3), and (5), it suffices to find either $\mathbf{f}$ or $\mathbf{h}$. To see this, note that if $\mathbf{f}$ is known, then $\mathbf{h}$ can be calculated from (3)–(5) and $h_{r}$. On the other hand, if $\mathbf{h}$ is known, the flows $\mathbf{f}$ can be found thanks to the monotonicity of (3) and the monotonicity of (5) in $[\underline{f}_{mn},\bar{f}_{mn}]$. In a nutshell, solving the WF task amounts to finding either $\mathbf{f}$ or $\mathbf{h}$. Because this simple observation is used throughout our analysis, it is summarized as a lemma.

Lemma 1**.**

Given $\mathbf{d}$, a triplet $(\mathbf{d},\mathbf{f},\mathbf{h})$ satisfying (2), (3), and (5) is uniquely characterized by $\mathbf{f}$ or $\mathbf{h}$.

The WF task can be posed as the feasibility problem

[TABLE]

Since (3)–(5) are quadratic equalities, problem (W1) is non-convex. For a tree WDS graph $\mathcal{G}$ (for which $P=N-1$), the edge-node incidence matrix $\mathbf{A}$ has $N-1$ linearly independent rows [24]. Then the flows $\mathbf{f}$ can be found uniquely from (2), and the WF task is readily solved according to Lemma 1. However, handling the WF task in a loopy $\mathcal{G}$ containing pumps remains non-trivial. Naturally, analyzing the uniqueness of a WF solution is a critical task. Reference [11] establishes the uniqueness of a WF solution under the assumption that the pressure gain added by a pump is constant. The next claim extends this uniqueness result regardless of the structure of $\mathcal{G}$ and to the more detailed pump model of  (5).

Theorem 1**.**

If the WF equations are feasible for some $\mathbf{d}$, they feature a unique solution.

Proof:

Proving by contradiction, assume $(\mathbf{h},\mathbf{f})$ and $(\tilde{\mathbf{h}},\tilde{\mathbf{f}})$ are two distinct solutions of (W1). Since the two flow vectors $\mathbf{f}$ and $\tilde{\mathbf{f}}$ satisfy (2), the vector $\mathbf{n}:=\tilde{\mathbf{f}}-\mathbf{f}$ lies in the nullspace of $\mathbf{A}^{\top}$, that is $\mathbf{A}^{\top}\mathbf{n}=\mathbf{0}$. Then it follows that

[TABLE]

Let us decompose $\mathbf{n}$ into its positive and negative entries as $\mathbf{n}=\mathbf{n}_{+}-\mathbf{n}_{-}$, where $\mathbf{n}_{+}\geq\mathbf{0}$ and $\mathbf{n}_{-}\geq\mathbf{0}$. Plugging this decomposition into (6) yields

[TABLE]

Recall from (3) that the pressure drop across a lossy pipe is monotonically increasing in flow. Similarly, the pressure added by a pump described by (5) is monotonically decreasing in flow. In other words, the pressure drop along a pump is monotonically increasing in flow. Hence, for each edge $p\in\mathcal{P}$:

[TABLE]

where $\mathbf{a}_{p}^{\top}$ is the $p$-th row of matrix $\mathbf{A}$.

Consider the $p$-th entry of vector $\mathbf{n}$. If $\tilde{f}_{p}>f_{p}$, then $n_{+,p}>0$ and $n_{-,p}=0$ by definition of $\mathbf{n}_{+}$ and $\mathbf{n}_{-}$. Then edge $p$ contributes to the left-hand side (LHS) of (7). The monotonicity in (8) entails that $\mathbf{a}_{p}^{\top}(\tilde{\mathbf{h}}-\mathbf{h})>0$ and so $n_{+,p}\cdot\mathbf{a}_{p}^{\top}(\tilde{\mathbf{h}}-\mathbf{h})>0$. The latter holds for all edges contributing to the LHS of (8), so that

[TABLE]

On the other hand, if $\tilde{f}_{p}<f_{p}$, then $n_{+,p}=0$ and $n_{-,p}>0$. Then, edge $p$ contributes to the right-hand side (RHS) of (8). The monotonicity in (8) entails $n_{-,p}\cdot\mathbf{a}_{p}^{\top}(\tilde{\mathbf{h}}-\mathbf{h})<0$. Applying the claim over all edges participating in the RHS of (8) provides

[TABLE]

The signs of the LHS and RHS contradict the equality in (7), thus proving the claim. ∎

Regarding the pressure drop law of (3), the Hazen-Williams equation is sometimes used wherein the flow $f_{mn}$ is raised to the exponent of $1.852$. This exponent is different from the exponent of $2$ in the Darcy-Weisbach equation; see for example [15]. While the Darcy-Weisbach equation is a theoretical formula, the Hazen-Williams equation is based on curve fitting of experimental data [4]. Either way, the uniqueness argument of Theorem 1 holds for any positive exponent on the water flow $f_{mn}$ involved in the pressure drop equation of (3).

The WF problem posed in Definition 1 follows the most traditional setup pursued in WF literature [11], [25], [26]. However, with multiple sources feeding a WDS an alternate form of WF becomes pertinent. Specifically, multiple supply nodes may operate as fixed-pressure nodes while the remaining nodes act as fixed-injection nodes [27]. The ensuing result extends the uniqueness of a WF solution to this more general setting.

Proposition 1**.**

Given: i) pressures $h_{i}$ for $i$ in a subset of nodes $\mathcal{N}_{r}\subset\mathcal{N}$; ii) injections $d_{i}$ for all $i\in\mathcal{N}\setminus\mathcal{N}_{r}$; and iii) the statuses and speeds for all pumps; the set of equations (2), (3), and (5) has at most one solution.

Proof:

Consider two pairs $(\mathbf{d},\mathbf{f})$ and $(\mathbf{d}^{\prime},\mathbf{f}^{\prime})$ satisfying (2). If $\mathbf{d}=\mathbf{d}^{\prime}$, the uniqueness of the WF solution follows from Theorem 1. Otherwise, [28, Lemma 3] dictates that there exists a path $\mathcal{P}_{mn}$ between fixed-pressure nodes $m$ and $n$, along which the flows in $\mathbf{f}$ differ from those in $\mathbf{f}^{\prime}$ in a consistent direction, i.e., either $f_{p}>f_{p}^{\prime}$ or $f_{p}<f_{p}^{\prime}$ for all edges $p\in\mathcal{P}_{mn}$. However, in that case, the monotonicity argument used in (8) would entail for the respective pressures $h_{m}-h_{n}\neq h_{m}^{\prime}-h_{n}^{\prime}$, which contradicts the hypothesis of $m$ and $n$ being fixed-pressure nodes. ∎

Having established the uniqueness of the WF solution, the ensuing sections develop a suite of WF solvers of complementary value: Each solver finds provably the WF solution under different network setups. The devised solvers exhibit different complexity, yet none of them requires a proper initialization. Figure 1 summarizes the particular WDS setups each solver can handle, and serves as a roadmap for the following sections. The solvers developed in this work are for the WF task as posed in Definition 1, whereas WF solvers for the setting considered in Proposition 1 are beyond the scope of this work.

IV Energy Function-Based Water Flow Solvers

For a WDS without pumps, the WF task simplifies to solving (2)–(3). These equations are structurally similar to the equations governing the flows and pressures in a natural gas network under steady-state conditions [28], [29]. Reference [25] finds the solution to the nonlinear WF equations in networks without pumps as the minimizers of a convex energy function; see also [26] and [30] for counterparts in gas networks. Upon briefly reviewing these reformulations, this section expands their scope to WDS with no pumps in cycles and other pressure drop laws; see roadmap of Figure 1.

IV-A Finding the WF Solution in WDS without Pumps

In the absence of pumps, the WF problem (W1) has been posed as a constrained minimization over $\mathbf{f}$ [25].

Lemma 2** ([25]).**

In a WDS without pumps ($\mathcal{P}=\bar{\mathcal{P}}_{a}$), the vector of water flows $\mathbf{f}$ satisfying (2)–(3) can be found as the unique minimizer of the convex minimization problem

[TABLE]

Moreover, if $\boldsymbol{\xi}^{*}$ is the vector of optimal Lagrange multipliers corresponding to (9b), then the nodal pressures are provided by $\mathbf{h}=(\mathbf{I}-\mathbf{1}\mathbf{e}_{r}^{\top})\boldsymbol{\xi}^{*}+h_{r}\mathbf{1}$.

The claim follows from the optimality conditions of (9). The strict convexity of (9a) shows that the WF solution is unique in networks without pumps. This claim has also been established using a contraction argument in [15]. Theorem 1 and Proposition 1 generalize this uniqueness claim to networks having pumps and multiple fixed-pressure nodes.

We next present two ways for handling (9): an second-order cone program (SOCP) reformulation, and a dual decomposition approach. It is known that $\ell_{p}$-norms can be handled using a hierarchy of second-order cone (SOC) constraints; see [31, Eq. (11)]. In the context of WDS, a more compact hierarchy of SOCs was originally put forth in [6] to handle the cubic powers of (9a) involved in an OWF problem. Adopting the latter reformulation, the minimizer of (9) can be found by solving the SOCP

[TABLE]

The fact that the $f_{p}$’s minimizing (9) and (10) coincide can be established by showing that the solution of (10) satisfies $w_{p}=|f_{p}|$; $y_{p}=w_{p}^{2}=f_{p}^{2}$; and $t_{p}=w_{p}^{3}=|f_{p}|^{3}$ for all $p\in\mathcal{P}$; see also [6]–[7] for details. Each constraint in (10d) is a rotated second-order cone [31].

Alternatively, problem (9) can be tackled via dual decomposition: During its $k$-th iteration, the primal variable $\mathbf{f}$ is updated by minimizing the Lagrangian function $L(\mathbf{f};\boldsymbol{\xi}^{k})$ evaluated at the latest estimate of dual variables $\boldsymbol{\xi}^{k}$. The latter problem decouples across pipes and enjoys a closed-form solution as

[TABLE]

For a step-size $\mu>0$, the dual variables are updated as

[TABLE]

Again for WDS without pumps, the WF task can alternatively be posed as an unconstrained minimization over $\mathbf{h}$ [25]

[TABLE]

The objective of (11) is convex (by composition rules) and differentiable. The $n$-th entry of its gradient $\mathbf{g}_{u}(\mathbf{h})$ is

[TABLE]

Setting this gradient equal to zero yields the WF equations after eliminating $\mathbf{f}$ from (2)–(3). The ambiguity in pressures can be waived by shifting the minimizer of (11) to match the reference pressure $h_{r}$. Once the pressures $\mathbf{h}$ are found, the flow vector $\mathbf{f}$ can be retrieved by Lemma 1. Problem (11) is amenable to any first-order method for unconstrained optimization, such as the gradient descent iterations $\mathbf{h}^{k+1}:=\mathbf{h}^{k}-\mu\mathbf{g}_{u}(\mathbf{h}^{k})$ for a step-size $\mu>0$, or accelerated variants.

IV-B Extension to WDS with no Pumps in Cycles

Since problems (9) and (11) cannot handle water networks with pumps, this section extends their applicability to networks with pumps, but not on cycles. This is accomplished by adopting the reduction technique of [32] to build what we term stitching algorithm:

1. S1)

Remove the edges of $\mathcal{G}$ corresponding to pumps $\mathcal{P}_{a}$. The obtained graph contains $|\mathcal{P}_{a}|+1$ disconnected components $\mathcal{G}_{c}$ for $c=1,\ldots,|\mathcal{P}_{a}|+1$. 2. S2)

Replace each component $\mathcal{G}_{c}$ by a supernode, and connect the supernodes using the edges in $\mathcal{P}_{a}$ to create the supergraph $\mathcal{G}^{\prime}$. Graph $\mathcal{G}^{\prime}$ features a tree structure. 3. S3)

Find the total water injection per supernode $\mathcal{G}_{c}$. Since $\mathcal{G}^{\prime}$ is a tree, the water flows on $\mathcal{P}_{a}$ can be found readily. 4. S4)

If the flow along pump $(m,n)\in\mathcal{P}_{a}$ is $f_{mn}$, modify the injections at nodes $m$ and $n$ as $\hat{d}_{m}=d_{m}-f_{mn}$ and $\hat{d}_{n}=d_{n}+f_{mn}$. 5. S5)

Solve (9) per connected component $\mathcal{G}_{c}$ to find the water flows at all lossy pipes. 6. S6)

Having acquired vector $\mathbf{f}$, the pressure vector $\mathbf{h}$ can be found from Lemma 1.

In S5), rather than solving the constrained minimization of (9), one could use its unconstrained counterpart of (11). Steps S1)–S4) are still needed to find a meaningful water injection vector per component $\mathcal{G}_{c}$. Once a pressure vector has been found per component, the pressures must be revised as follows: The pressures within the component containing the reference node, say component $\mathcal{G}_{1}$, are kept unaltered. Consider a pump running from node $m\in\mathcal{G}_{1}$ to node $n\in\mathcal{G}_{2}$. Using the flow $f_{mn}$ computed in step S3), the pressure gain $g_{mn}=g_{mn}(f_{mn};\omega^{0}_{mn})$ can be found from (5). Having found $h_{m}$ by solving (11) in $\mathcal{G}_{1}$, the pressure at node $n$ can be computed as $h_{n}=h_{m}+g_{mn}$. If the pressure at node $n$ recovered by solving (11) in $\mathcal{G}_{2}$ is $\tilde{h}_{n}$, then all pressures within $\mathcal{G}_{2}$ should be shifted by $h_{n}-\tilde{h}_{n}$. The process is repeated for all components to recover the entire pressure vector $\mathbf{h}$.

IV-C Handling the Hazen-Williams Pressure Drop Law

The aforementioned WF solvers can be modified to handle the Hazen-Williams in lieu of the Darcy-Weisbach pressure drop equations. Recall that according to the Hazen-Williams pressure drop law, the exponent of $|f_{mn}|$ in (3) changes from $2$ to $\rho=1.852$. To handle this, problem (9) is generalized to

[TABLE]

where the exponent $\rho$ depends on the pressure drop law and could vary per pipe. Problem (12) remains convex for $\rho>0$.

Even though the SOCP model of (10) cannot be used to solve (12), a different SOCP reformulation can be derived as long as $\rho$ is rational; see [31, pp. 12–13] for details. Moreover, the dual decomposition updates of Section IV-A can be adapted to solve (12). Similarly, the unconstrained minimization of (11) can be modified to

[TABLE]

Heed that the energy function-based solvers for (W1) can not handle the WF setups with multiple fixed-pressure nodes described under Proposition 1. This is because energy function-based solvers do not take reference pressure as an input to the problem. Rather, the obtained pressures are shifted a posteriori to agree with a given fixed pressure. For a single fixed-pressure node, the applicability of our WF solvers can be broadened to WDS with pumps in cycles. To do so, we next pursue an MI-QCQP solver. Unlike the energy function-based WF solvers of this section, the MI-QCQP solver of Section V applies only to the Darcy-Weisbach equation; its possible extension to the Hazen-Williams pressure drop law goes beyond the scope of this work.

V MI-QCQP Water Flow Solver

This section presents an MI-QCQP relaxation of (W1) along with sufficient conditions for its exactness.

V-A Problem Reformulation

The non-convexity of the WF problem (W1) is due to the quadratic equality constraints of (3) and (5). Adopting the convex relaxation approach of [7]–[6], the pressure drop along pipe $(m,n)\in\bar{\mathcal{P}}_{a}$ is relaxed from (3) to

• •

$h_{m}-h_{n}\geq c_{mn}f_{mn}^{2}$ for $f_{mn}\geq 0$; or

• •

$h_{n}-h_{m}\geq c_{mn}f_{mn}^{2}$ for $f_{mn}\leq 0$.

To handle the two cases, let us introduce a binary variable $x_{mn}$ capturing the flow direction on pipe $(m,n)$. By using the so-termed big-$M$ trick, the two cases can be modeled as

[TABLE]

The same MI-QCQP model has been also advocated for handling the pressure drop equations along the pipelines in natural gas networks; see e.g., [33], [34], [28]. Constraint (14a) implies that $x_{mn}=\operatorname{sign}(f_{mn})$. If $x_{mn}=1$, the convex quadratic constraint in (14b) is activated, whereas constraint (14c) becomes trivial. The claim reverses for $x_{mn}=0$. Depending on the value $x_{mn}$, if either (14b) or (14c) is satisfied with equality, the relaxation is deemed as exact. If that happens for all lossy pipes, the feasible set of (14) has captured the original non-convex constraints in (3).

Similarly, the pressure added by pump $(m,n)\in\mathcal{P}_{a}$ is relaxed from (5) to

[TABLE]

which is a convex constraint since $\lambda_{mn}<0$. Again, if the inequality in (15) is satisfied with equality for all pumps, the relaxation of (5) to (15) is deemed as exact.

One may now try solving the feasibility problem in (W1) by replacing (3)–(5) by (14)–(15). If all the relaxed constraints are satisfied with equality, the relaxation is exact. Unfortunately, toy numerical tests can verify that such relaxation is inexact even for simple water networks. To favor exact relaxations, we convert the feasibility to a minimization problem with a judiciously selected cost. We pose the MI-QCQP

[TABLE]

where the binary vector $\mathbf{x}$ contains all $\{x_{mn}\}_{(m,n)\in\bar{\mathcal{P}}_{a}}$, and the cost is defined as

[TABLE]

The cost sums up the absolute pressure drops across lossy pipes minus the pressure gains added by pumps. This is in contrast to the MI-QCQP of [29] that recovers a gas flow solution, whose penalty did not include compressors. It will be shown in Section V-B that the MI-QCQP formulation (W2) guarantees exactness under realistic network conditions.

Problem (W2) is non-convex due to the binary variables $\mathbf{x}$. Given the advancements in MI-QCQP solvers, this minimization can be handled for moderately sized networks. Recall that every lossy pipe is associated with one binary and one continuous optimization variable. To accelerate the computations, we next provide a simple pre-processing step to determine the flows in all edges not belonging to a cycle.

Lemma 3**.**

Let $\mathbf{f}$ be the unique flow solution of (W1), and consider any vector $\tilde{\mathbf{f}}$ satisfying $\mathbf{A}^{\top}\tilde{\mathbf{f}}=\mathbf{d}$. For any edge $(m,n)$ not belonging to a cycle, it holds that $\tilde{f}_{mn}=f_{mn}$.

Proof:

Consider the minimum-norm solution $\mathbf{f}_{0}:=\left(\mathbf{A}^{\top}\right)^{\dagger}\mathbf{d}$ to the linear system of (2), where $\left(\mathbf{A}^{\top}\right)^{\dagger}$ is the pseudo-inverse of $\mathbf{A}^{\top}$. Any other solution of (2) can be expressed as $\mathbf{f}=\mathbf{f}_{0}+\mathbf{n}$ for some vector $\mathbf{n}\in\operatorname{null}(\mathbf{A}^{\top})$.

The space $\operatorname{null}(\mathbf{A}^{\top})$ can be represented using the cycles of graph $\mathcal{G}$ as explained next. For cycle $\mathcal{C}$ in $\mathcal{G}$, select an arbitrary direction and define its indicator vector $\mathbf{n}_{\mathcal{C}}\in\mathbb{R}^{P}$ as

[TABLE]

We are particularly interested in a set of fundamental cycles defined as follows [24]: In graph $\mathcal{G}=(\mathcal{N},\mathcal{P})$, select a spanning tree $\mathcal{T}$. Every edge $p\in\mathcal{P}\setminus\mathcal{T}$ along with some edges of $\mathcal{T}$ form a cycle. The cycles formed using all edges in $\mathcal{P}\setminus\mathcal{T}$ comprise the set $\mathcal{L}_{\mathcal{T}}$ of fundamental cycles. In general, a graph may have multiple spanning trees, making the set of fundamental cycles non-unique.

A key property of fundamental cycles is that the space $\operatorname{null}(\mathbf{A}^{\top})$ is spanned by the indicator vectors for any set of fundamental cycles [24, Corollary 14.2.3]. Therefore, the vector $\mathbf{n}$ in $\mathbf{f}=\mathbf{f}_{0}+\mathbf{n}$ can be decomposed as

[TABLE]

where $\mathbf{n}_{\ell}$ is the indicator vector for cycle $\ell$ and $\alpha_{\ell}\in\mathbb{R}~{}\forall\ell$.

Consider the $p$-th entry of $\mathbf{n}$. If edge $p\in\mathcal{P}$ does not belong to any cycle, then it does not belong to any fundamental cycle. Then $n_{\ell,p}=0$ for all $\ell\in\mathcal{L}_{\mathcal{T}}$, and $n_{p}=0$ follows from (16). To conclude, if edge $p$ does not belong to a cycle, then all solutions $\mathbf{f}$ of (2) agree in their $p$-th entry. ∎

Lemma 3 ensures that for weakly meshed WDS, the number of variables in (W2) can be reduced markedly. Moreover, the flow on any edge not belonging to a cycle coincides with the related entry of the minimum-norm solution to (2).

V-B Exactness of the Relaxation

The next result provides conditions under which a minimizer of (W2) satisfies (14)–(15) with equality.

Theorem 2**.**

In a WDS where no edge $(m,n)\in\mathcal{P}$ belongs to more than one cycle, a minimizer of (W2) minimizes (W1) as well, if (W1) is feasible.

Theorem 2 (shown in the appendix) asserts that the MI-QCQP relaxation of (W1) to (W2) is exact in water networks with non-overlapping cycles, that is cycles sharing no edges. The claim holds regardless of water demand or the presence of pumps. This is in contrast to the exactness claims of [29], where compressors were not allowed on cycles. Although the assumptions of Theorem 2 may not always hold, the numerical tests of Section VII indicate the relaxation is exact and (W2) solves (W1) for practical WDS with overlapping cycles.

The condition of Theorem 2 cannot be relaxed analytically: One could construct pathological WDS with overlapping cycles that render the relaxation of (W2) inexact. To present such a counterexample, consider the $4$-node and $5$-pipe WDS of Fig. 2. The coefficients $c_{mn}$’s are shown on the respective pipes. Nodes $2$ and $4$ host water demands, and the reference node $1$ supplies water at the reference pressure of $h_{1}=10$. The edge $(1,3)$ belongs to two cycles and hence the condition of Theorem 2 is violated. The minimizer of (W2) satisfies the relaxed constraint (14) with strict inequality. The pressures at nodes $1$ and $3$ were found to differ by $h_{1}-h_{3}=3.435$, while the frictional drop was $c_{13}f_{13}^{2}=0.014$. This WDS though is characterized by a strong disparity of pipe coefficients. For a WDS with limited variations in pipe dimensions and material, such disparity is not anticipated in practice. Similar to the example of Fig. 2, one can construct setups with multiple fixed-pressure nodes for which the MI-QCQP relaxation is inexact.

Since the WDS of Fig. 2 does not host pumps, the solution to the previous WF problem was eventually found via (9). This motivates us to exploit the ability of the energy function-based approach to handle overlapping cycles alongside the merit of the MI-QCQP relaxation to handle pumps in cycles. To solve the WF problem for a broader class of WDS network topologies, a hybrid solver is deferred to Section VI. Before that, the formulation of (W2) is contrasted to prior work.

V-C Comparison to Prior Work and Extensions

The proposed WF solver involves three ingredients:

1. I1)

replacing (3) by the MIQP model of (14); 2. I2)

replacing (4) by the convex constraint of (15); and 3. I3)

solving (W2) in lieu of (W1).

The key contribution of this section is I3), since I1) and I2) have appeared before: References [6]–[7] and [35] consider the OWF problem assuming known flow directions (say $f_{mn}\geq 0$ for all pipes), and relax (3) to [cf. (14)]

[TABLE]

Regarding pumps, references [6]–[7] and [35] consider only variable-speed pumps. Moreover, each pump speed $\omega_{mn}$ is limited within $0\leq\omega_{mn}\leq\bar{\omega}_{mn}$, which may be unrealistic since pumps come with positive lower speed bounds. Substituting this speed range into the monotonic formula of (5) yields the convex quadratic constraint

[TABLE]

Heed that the actual speed $\omega_{mn}$ has been eliminated from (18). Nevertheless, once an OWF minimizer is found [and so a triplet $(h_{m},h_{n},f_{mn})$ satisfying (18)], the speed $\omega_{mn}$ making (18) an equality can be readily recovered [6].

References [6]–[7] aim at minimizing the energy losses across pipes. Ignoring the details of geographical elevation and electricity prices, their cost simplifies to $\sum_{(m,n)\in\bar{\mathcal{P}}_{a}}c_{mn}f_{mn}^{3}$. Thanks to its form, minimizing this cost subject to (17)–(18) renders the relaxation in (17) numerically exact. This was also analytically shown granted the simplifying assumptions of: a1) no pumps on cycles; a2) known flow directions; a3) if a node is fed by more than one pipes, all of them have to be equipped with valves; and a4) no fixed-speed pumps. Reference [35] considers the joint scheduling of WDS and electric power distribution networks. Since the exactness of (17) is not promoted by the cost anymore, a feasible point pursuit-based approach is applied to find a stationary point satisfying (17) with equality.

From the preceding review, it is evident that albeit a relaxation of (3) is convenient, its exactness is not always guaranteed. While some OWF instances feature exactness [6]–[7]; others call for additional measures [35]. We shall next address the natural question: Q1) Can one use one of the aforesaid OWF solvers to solve the WF task? Numerical tests on WF instances even in simple WDS demonstrate that the approaches of [6]–[7] and [35] fail to solve (W1), that is the convex relaxation is not exact. This is because even if an instance of (W1) satisfies all other conditions needed for each one of these solvers, the common assumption a4) is violated since pump speeds are fixed for (W1).

One may reverse question Q1) and pose Q2): Can the developed WF solvers be used towards solving OWF tasks? One could identify two possible uses. First, the proposed tools can be used instead of EPANET in existing zero-order OWF algorithms that rely on a WF solver. Second, the success of the proposed penalized relaxation could be adopted in OWF tasks. For example, our previous work [8] deals with optimal pump scheduling under dynamic electricity pricing. This OWF task handles unknown flow directions and is relaxed to an MI-SOCP after $g_{mn}(f_{mn};\omega^{0}_{mn})$ in (5) is approximated as constant. The relaxation is provably exact under certain conditions and after adding a penalty to the cost of the MI-SOCP. The relaxation is numerically exact under a broader range of WDS conditions. Interestingly, the penalty that worked for [8] does not work for the WF task studied here. For this reason, this work put forth the penalty $s(\mathbf{h})$ in (W2) and took a totally different route for establishing exactness. The sufficient conditions for exactness of Theorem 2 are significantly simpler compared to those in prior works on OWF.

VI A Hybrid Water Flow Solver

The hybrid WF solver relaxes the assumption of non-overlapping cycles of Theorem 2 to the following condition.

Assumption 1**.**

The water network has no pumps in overlapping cycles.

The assumption permits overlapping cycles, but these particular cycles should carry no pumps. For a network satisfying Assumption 1, the WF task can be solved through the following steps illustrated also in Fig. 3:

1. T1)

Any cycle $\mathcal{C}$ with pumps is non-overlapping. Thus, for any node $n$ belonging to $\mathcal{C}$, two cases arise:

i) Node $n$ belongs only to cycle $\mathcal{C}$ (node $3$ of $\mathcal{C}_{1}$); or

ii) Node $n$ belongs to other cycle(s) $\mathcal{C}^{\prime}$ as well; yet $\mathcal{C}$ and $\mathcal{C}^{\prime}$ share no edge (e.g., node $4$ of $\mathcal{C}_{2}$).

Then, reduce graph $\mathcal{G}$ to $\mathcal{G}^{\prime}$ through the steps:

• •

If for a cycle $\mathcal{C}$ having pumps, all nodes do not belong to any other cycle, replace $\mathcal{C}$ by a supernode $n_{\mathcal{C}}$.

• •

If for a cycle $\mathcal{C}$ having pumps, node $n$ belongs to other cycle(s), identify the edges $(n_{1},n)$ and $(n,n_{2})$, which belong to $\mathcal{C}$. These two edges belong to no other cycles as $\mathcal{C}$ is non-overlapping. Split the node $n$ into $n$ and $n^{\prime}$ connected by a lossless edge $(n,n^{\prime})$, such that all edges in $\mathcal{G}$ other than $(n_{1},n)$ and $(n,n_{2})$ that were incident on $n$ are now incident on $n^{\prime}$. After repeating this step for all nodes in $\mathcal{C}$ that belong to multiple cycles, replace $\mathcal{C}$ by a supernode $n_{\mathcal{C}}$.

This process ensures that in $\mathcal{G}^{\prime}$ all supernodes and the pumps left out from $\mathcal{G}$ do not appear on cycles. 2. T2)

Use Lemma 3 to compute the water flows on the edges $(m,n_{\mathcal{C}})$ incident to supernodes $n_{\mathcal{C}}$’s in $\mathcal{G}^{\prime}$. 3. T3)

For each edge $(m,n_{\mathcal{C}})$ in $\mathcal{G}^{\prime}$, modify the injection at node $m$ as $\hat{d}_{m}=d_{m}-f_{mn_{\mathcal{C}}}$. 4. T4)

Partition $\mathcal{G}^{\prime}$ into a set of connected components $\mathcal{G}_{c}$ by removing the supernodes and their incident edges. Each $\mathcal{G}_{c}$ has known injections and bears no pumps in cycles. 5. T5)

Use the stitching algorithm of Section IV-B per component $\mathcal{G}_{c}$ to find the flows within $\mathcal{G}_{c}$. 6. T6)

Each supernode $n_{\mathcal{C}}$ is split back to the cycle $\mathcal{C}$ it replaced. If the edge $(m,n_{\mathcal{C}})$ of $\mathcal{G}^{\prime}$ corresponded to edge $(m,n)$ of $\mathcal{G}$, modify the injection at node $n\in\mathcal{C}$ as $\hat{d}_{n}=d_{n}+f_{mn_{\mathcal{C}}}$. 7. T7)

Solve (W2) per cycle $\mathcal{C}$ to find the water flows on $\mathcal{C}$. 8. T8)

Given vector $\mathbf{f}$, find the pressure vector $\mathbf{h}$ using Lemma 1.

Let us apply the previous steps on the $23$-node and $28$-pipe water network of Fig. 3. There are two cycles carrying pumps, marked as $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$. Node $4$ of cycle $\mathcal{C}_{1}$ belongs to multiple cycles and hence it is split in $4$ and $4^{\prime}$. Next, removing $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ along with the edges that connect these cycles with the rest of the graph results in the connected subgraphs $\mathcal{G}_{1}$, $\mathcal{G}_{2}$, and $\mathcal{G}_{3}$ shown on the bottom of Fig. 3. The demands on boundary nodes are modified as per steps T3) and T6). The edge flows within each connected component are subsequently found by solving (9). The flows on $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ are finally found by (W2).

This WDS setup could not be handled by the energy function-based approach of Section IV alone due to the presence of pumps on cycles $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$. The convex relaxation of Section V alone is not guaranteed to succeed either, due to the presence of overlapping cycles. Combining the merits of each method and leveraging Lemma 3, this hybrid method can handle successfully this WDS setup.

VII Numerical Tests

The new WF solvers were evaluated under the different network configurations of Fig. 1. First, the new WF solvers were evaluated on the EPANET Example Network-2 representing a water network from Cherry Hills, Connecticut shown in Fig. 4 [2]. This network consists of $P=40$ pipes, $N=34$ demand nodes, one tank, and one pump station. We modified the network by representing the pump station as a reservoir with pressure $100$ ft connected to a fixed-speed pump. All nodes were assumed at the same elevation. The pipe friction coefficients $c_{mn}$’s, and base demand vector $\mathbf{d}$, were derived from the EPANET benchmark. Note that the WDS in Fig. 4 has overlapping cycles but the pump is not in a loop. Thus, it qualifies for being solved using the energy function method alongside the stitching algorithm of Section IV-B.

We tested whether the proposed solvers yield the same solution as EPANET. For this purpose, we obtained WF solutions using the constrained energy function minimization of (9) and the MI-QCQP of (W2). Problem (9) was solved using the closed-form dual decomposition steps of Section IV-A for $\mu=10^{-4}$ within $20,000$ iterations. The MI-QCQP in (W2) was solved using the MATLAB-based toolbox YALMIP along with the mixed-integer solver CPLEX [36], [37]. All tests were run on a $2.7$ GHz, Intel Core i5 computer with $8$ GB RAM. The flows obtained by the two solvers were very close to the EPANET solution, as illustrated in Fig. 5. EPANET uses more detailed flow models, e.g., the coefficient $c_{mn}$ in (3) depends weekly on flow $f_{mn}$. The $20,000$ iterations for the primal-dual updates of (9) were completed within $4.5$ sec, and the running time of the MI-QCQP (W2) was $1.6$ sec.

We next evaluated the performance of the MI-QCQP-based solver of (W2) in finding the correct WF solution, when the conditions of Theorem 2 are not satisfied. Specifically, the network of Fig. 4 has overlapping cycles. To represent various demand levels, we generated $100$ random WF instances by scaling the benchmark demand $\mathbf{d}$ by a scalar uniformly drawn from $[0,1.5]$. Given a minimizer of (W2), we defined the inexactness over lossy pipe $(m,n)$ as ${|h_{m}-h_{n}|-c_{mn}f_{mn}^{2}}$. For each run of (W2) with a random input, the ranked maximum inexactness gap over all lossy pipes is displayed on Fig. 6 (top). Despite violating the conditions of Th. 2, the inexactness gap was small for all random tests. Computationally, the running time of (W2) over the $100$ instances with a median value of only $1.68$ sec; see Fig. 6 (bottom).

Our WF solvers were subsequently tested on the 20-node synthetic network of Fig. 7 and under different demands. The WDS of Fig. 7 was created by combining two copies of the popular 10-node WDS representing the modified Arava Valley network from Israel; see [8], [22] for network data. The additional pipes $(1,2)$, $(7,14)$, and $(8,20)$ have the same dimensions as pipes $(4,5)$, $(5,6)$, and $(6,8)$, respectively. These were added to generate different network conditions. The parameters for all the pumps were derived by scaling the pump parameters of [7] by $0.25$, and were set to $(\lambda_{mn},\bar{\mu}_{mn},\bar{\nu}_{mn})=(-2.735\cdot 10^{-5},0.0129,55.83)$.

The performance of (W2) was tested for two network configurations derived from the WDS of Fig. 7: C-i) pipe $(4,10)$ removed; C-ii) pipe $(8,20)$ removed. Note that configuration C-i) has pumps in non-overlapping cycle along side other overlapping cycles not hosting pumps. Thus, the WF problem for C-i) can provably be solved using the hybrid WF solver of Section VI. Moreover, configuration C-ii) has pumps in overlapping cycles and thus cannot be provably solved even by the hybrid WF solver. However, we numerically tried (W2) in solving $100$ random WF instances for both C-i) and C-ii).

The $100$ random WF instances were generated by fixing pressure $h_{1}=10$ m and drawing the remaining pressures independently as Gaussian random variables of mean $10$ m and variance $2$ m2, or $\mathcal{N}(10,2)$. The on/off statuses of pumps were drawn as independent Bernoulli random variables with mean of $0.5$. The flow for on pumps was drawn with uniform probability on the allowable pump flow $[250,1500]$ m${}^{3}/$hr. For off pumps, flows were drawn from $\mathcal{N}(0,200)$. The pressure added by pumps was then calculated using (5). For pump $(m,n)$, the receiving node pressure was updated as $h_{n}:=h_{m}+g_{mn}$. The water flow in all lossy pipes was calculated from the so obtained pressures and (3). Once the complete flow vector was obtained, the injections were computed from (21). The obtained vector of injections, pump on/off status, and pressure $h_{1}$ served as a feasible input for the WF task. Similarly $100$ feasible random WF instances were generated for configuration C-ii), and solved using (W2). The maximum time for solving (W2) was set to $1$ min in CPLEX.

The WF task for both configurations was solved within a $1$-minute deadline for $98$ out of the $100$ random instances with the average time being $1.01$ and $1.06$ sec for C-i) and C-ii), respectively. The ranked maximum inexactness gap over all lossy pipes for C-i) and C-ii) is shown in Fig. 8. Although the conditions of Th. 2 were not satisfied by either configuration, the maximum inexactness obtained was less than $10^{-3}$ for $74$ and $92$ out of the $98$ instances of C-i) and C-ii that were solved. For the two cases where the MI-QCQP failed to converge within the $1$-min deadline, we modified the big-$M$ parameter value from $M=300$ to $M=80$. The result was that the two instances were solved in $0.13$ and $0.82$ seconds with the inexactness gap being less than $10^{-5}$.

To evaluate whether the status of the pumps has any significant effect on runtime, we conducted the previous tests on the WDS of Fig. 7 with all the pumps running. The inexactness gap and running times were found similar to the tests with pump statuses randomly chosen. Interestingly, the inexactness gap and running times for the new tests did not show significant differences from the tests with random pump statuses; see Fig. 9. 99 out of the $100$ random WF instances were solved within the $1$-min deadline with a median running time of $0.75$ sec. Thus, the MI-QCQP (W2) is in general a powerful WF solver for various network configurations.

VIII Conclusions

Using recent tools from graph theory, convex relaxations, energy function-based approaches, and mixed-integer programing, this work has provided a fresh perspective on the physical laws governing water distribution networks. It has been established that the WF problem admits a unique solution even in networks with multiple fixed-pressure nodes and flow-dependent pump models. This WF solution can be provably recovered via a hierarchical stack of WF solvers suitable for different network configurations. Radial networks can be handled by simple convex minimization tasks, whereas networks with cycles call for more elaborate MI-QCQP-based solvers. Nevertheless, numerical tests demonstrate that even the MI-QCQP approach scales well for moderately-sized networks. The MI-QCQP solvers have been derived upon a convex relaxation of the pressure drop equation followed by an objective penalization, an approach that sparked a parallel line of research on the OWF problem [8]. Network configurations hosting pumps on overlapping cycles remain to be a challenging case.

Proof:

The cost of (W2) can be written as

[TABLE]

To express pressure differences along the flow direction in a compact manner, define the $P\times N$ edge-node incidence matrix $\mathbf{A}(\mathbf{f})$: Its dependence on $\mathbf{f}$ signifies that the directionality of each edge coincides with the flow directions in $\mathbf{f}$. Therefore, if the $p$-th row of $\mathbf{A}(\mathbf{f})$ is associated with pipe $p=(m,n)$, then its $(p,k)$ entry is

[TABLE]

For zero flows $(\operatorname{sign}(f_{mn})=0)$, the default pipe direction $(m,n)$ is selected without loss of generality.

Based on $\mathbf{A}(\mathbf{f})$, the pressure differences along the direction of flows can be written as $\mathbf{A}(\mathbf{f})\mathbf{h}$, and so

[TABLE]

If (W1) is feasible, denote its unique solution by $(\bar{\mathbf{f}},\bar{\mathbf{h}})$. Since the base directionality of the WDS graph $\mathcal{G}$ is arbitrary, let it coincide with the water flow directions of $\bar{\mathbf{f}}$. Using this convention, it follows that $\bar{\mathbf{f}}\geq\mathbf{0}$ and $\mathbf{A}(\bar{\mathbf{f}})$ is identical to the base incidence matrix $\mathbf{A}$. Next, proving by contradiction, suppose $(\tilde{\mathbf{f}},\tilde{\mathbf{h}})$ is a minimizer of (W2), which is not feasible for (W1). Since both $\bar{\mathbf{f}}$ and $\tilde{\mathbf{f}}$ satisfy (2) for $\mathbf{A}=\mathbf{A}(\bar{\mathbf{f}})$, there exists a nonzero vector $\mathbf{n}\in\operatorname{null}(\mathbf{A}^{\top})$ such that

[TABLE]

As shown in (16), vector $\mathbf{n}$ can be expressed in terms of a set $\mathcal{L}_{\mathcal{T}}$ of fundamental cycles using the flow directions of $\mathbf{A}$. To simplify notation, let us define $\tilde{\mathbf{A}}:=\mathbf{A}(\tilde{\mathbf{f}})$. Since every edge is assumed to belong to at most one cycle, the set of fundamental cycles is unique irrespective of the spanning tree $\mathcal{T}$ selected. Therefore, the set of fundamental cycles in fact contains all cycles in the graph, implying $\mathcal{L}_{\mathcal{T}}=\mathcal{L}$.

Building on (16), consider a decomposition of $\mathbf{n}$ as a weighted sum of indicator vectors $\mathbf{n}_{\ell}$’s. If there exists a cycle $\ell$, for which $\alpha_{\ell}<0$, then one can reverse the direction of cycle $\ell$ and substitute $(\alpha_{\ell},\mathbf{n}_{\ell})$ in (16) with $(-\alpha_{\ell},-\mathbf{n}_{\ell})$. Hence, it can be assumed that $\alpha_{\ell}\geq 0$ for all $\ell\in\mathcal{L}$ without loss of generality.

Each vector $\mathbf{n}_{\ell}$ can be decomposed as $\mathbf{n}_{\ell}=\mathbf{n}_{\ell}^{+}-\mathbf{n}_{\ell}^{-}$, where $\mathbf{n}_{\ell}^{+}:=\max\{\mathbf{n}_{\ell},\mathbf{0}\}$ and $\mathbf{n}_{\ell}^{-}:=\max\{-\mathbf{n}_{\ell},\mathbf{0}\}$. Because every edge belongs to at most one cycle, it holds that

[TABLE]

where the $p$-th entry of vector $\mathbf{n}^{0}$ is $1$ if edge $p$ does not belong to any cycle; and [math], otherwise. Using (21) in (19), the objective function of (W2) becomes

[TABLE]

Albeit $s(\mathbf{h};\mathbf{f})$ will be evaluated for different pairs $(\mathbf{h};\mathbf{f})$, the vectors $\mathbf{n}_{\ell}$’s remain unchanged and depend on $\mathbf{A}$. Based on (16) and (22), we will next show that $s(\bar{\mathbf{h}};\bar{\mathbf{f}})<s(\tilde{\mathbf{h}};\tilde{\mathbf{f}})$. To do so, we consider the three terms of (22) separately.

First summand of (22). Recall $\mathbf{n}_{\ell}^{+}$ is a binary vector; and the base graph directionality is such that $\bar{\mathbf{f}}\geq\mathbf{0}$. Consider the entries of $\bar{\mathbf{f}}$ and $\tilde{\mathbf{f}}$ for the edges $p$ related to $n_{\ell,p}^{+}=1$. If $n_{\ell,p}^{+}=1$, then $\tilde{f}_{p}=\bar{f}_{p}+\alpha_{\ell}\geq\bar{f}_{p}\geq 0$. In that case, if edge $p=(m,n)$ and relates to a lossy pipe, we get

[TABLE]

where the first inequality stems from constraint (14) of (W2); the second one from $\tilde{f}_{p}\geq\bar{f}_{p}\geq 0$; and the equality from constraint (3) of (W1).

If edge $p=(m,n)$ relates to a pump, we can also show

[TABLE]

Consider cycle $\ell$ and sum up the LHS and RHS of (23) or (VIII) for all $p$ with $n_{\ell,p}^{+}=1$ to get

[TABLE]

where the inequality is strict if $\alpha_{\ell}\neq 0$. Summing (25) over all cycles provides

[TABLE]

with strict inequality arising from the fact that not all $\alpha_{\ell}$’s can be zero for a nonzero $\mathbf{n}$.

Second summand of (22). As explained earlier, if $n_{\ell,p}^{+}=1$, then $\tilde{f}_{p}=\bar{f}_{p}+\alpha_{\ell}\geq\bar{f}_{p}\geq 0$ so $\tilde{f}_{p}$ remains positive. On the other hand, if $n_{\ell,p}^{-}=1$, then $\tilde{f}_{p}=\bar{f}_{p}-\alpha_{\ell}<\bar{f}_{p}$. In the latter case, the flow $\tilde{f}_{p}$ has decreased, and its sign may have been reversed to negative. Since all $\bar{f}_{p}$’s are positive, the flow reversals in $\tilde{f}_{p}$’s can be modeled as $\mathbf{A}=\mathbf{S}\tilde{\mathbf{A}}$, where matrix $\mathbf{S}:=\operatorname{dg}(\operatorname{sign}(\tilde{\mathbf{f}}))$ is diagonal with the signs of $\tilde{\mathbf{f}}$ on its main diagonal. Since the vectors $\mathbf{n}_{\ell}$’s form a basis for $\operatorname{null}(\mathbf{A}^{\top})$, it holds that $\mathbf{A}^{\top}\mathbf{n}_{\ell}=\mathbf{0}$ and so

[TABLE]

Similar properties hold for $\tilde{\mathbf{n}}_{\ell}:=\mathbf{S}\mathbf{n}_{\ell}$. To see this, the vector $\tilde{\mathbf{n}}_{\ell}$ belongs to $\operatorname{null}(\tilde{\mathbf{A}}^{\top})$ since

[TABLE]

Therefore, we also get that

[TABLE]

where $\tilde{\mathbf{n}}_{\ell}^{+}:=\max\{\tilde{\mathbf{n}}_{\ell},\mathbf{0}\}$ and $\tilde{\mathbf{n}}_{\ell}^{-}:=\max\{-\tilde{\mathbf{n}}_{\ell},\mathbf{0}\}$.

By definition of $\tilde{\mathbf{n}}_{\ell}$, if $n_{\ell,p}=1$, then $S_{p,p}=1$ and $\tilde{n}_{\ell,p}=1$. However, if $n_{\ell,p}=-1$, then $S_{p,p}=+1$ or $S_{p,p}=-1$ depending on the sign of $\tilde{f}_{p}$, and so $\tilde{n}_{\ell,p}=1$ or $\tilde{n}_{\ell,p}=-1$. It therefore follows that

[TABLE]

Back to the $\ell$-th term of the second summand in (22):

[TABLE]

where $(a)$ stems from (29b); $(b)$ from (28); $(c)$ from (29a); $(d)$ from (25); and $(e)$ from (27). Summing (VIII) over $\ell\in\mathcal{L}$

[TABLE]

where the strict inequality stems from that argument $(d)$ in (VIII) is strict for $\alpha_{\ell}\neq 0$, and not all $\alpha_{\ell}$’s can be zero.

Third summand of (22). The third summand sums up the pressure differences along the direction of flow for all edges not lying in any cycle. If edge $p=(m,n)$ belongs to this case (i.e., $n_{p}^{0}=1$), then $\bar{f}_{p}=\tilde{f}_{p}$ and so as in (23) we get

[TABLE]

Summing up over all edges with $n_{p}^{0}=1$ yields

[TABLE]

Adding (25), (31), and (32) by parts gives $s(\tilde{\mathbf{h}};\tilde{\mathbf{f}})>s(\bar{\mathbf{h}};\bar{\mathbf{f}})$. This contradicts that $(\tilde{\mathbf{f}},\tilde{\mathbf{h}})$ is a minimizer of (W2) since $(\bar{\mathbf{f}},\bar{\mathbf{h}})$ is feasible for (W2), and concludes the proof. ∎