Natural Gas Flow Solvers using Convex Relaxation

TL;DR

This paper introduces a convex relaxation-based solver for steady-state natural gas flow problems, proving its exactness under certain conditions and demonstrating its scalability and robustness through numerical tests.

Contribution

It develops a novel convex relaxation and a mixed-integer quadratic programming solver for natural gas flow equations, ensuring solution uniqueness and broad applicability.

Findings

The solution to the gas flow problem is unique under arbitrary network configurations.

The convex relaxation is exact in networks with non-overlapping cycles and a single reference node.

The MI-QCQP solver scales well and remains accurate even when sufficient conditions are not met.

Abstract

The vast infrastructure development, gas flow dynamics, and complex interdependence of gas with electric power networks call for advanced computational tools. Solving the equations relating gas injections to pressures and pipeline flows lies at the heart of natural gas network (NGN) operation, yet existing solvers require careful initialization and uniqueness has been an open question. In this context, this work considers the nonlinear steady-state version of the gas flow (GF) problem. It first establishes that the solution to the GF problem is unique under arbitrary NGN topologies, compressor types, and sets of specifications. For GF setups where pressure is specified on a single (reference) node and compressors do no appear in cycles, the GF task is posed as an convex minimization. To handle more general setups, a GF solver relying on a mixed-integer quadratically-constrained quadratic program (MI-QCQP) is also devised. This solver can be used for any GF setup at any NGN. It introduces binary variables to capture flow directions; relaxes the pressure drop equations to quadratic inequality constraints; and uses a carefully selected objective to promote the exactness of this relaxation. The relaxation is provably exact in NGNs with non-overlapping cycles and a single fixed-pressure node. The solver handles efficiently the involved bilinear terms through McCormick linearization. Numerical tests validate our claims, demonstrate that the MI-QCQP solver scales well, and that the relaxation is exact even when the sufficient conditions are violated, such as in NGNs with overlapping cycles and multiple fixed-pressure nodes.

Equations

$$q_n=\sum _{\ell :(n,k)\in \mathcal{P}}{\varphi }_{\ell }-\sum _{\ell :(k,n)\in \mathcal{P}}{\varphi }_{\ell }.$$

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

$${\mathbf{A}}^{\top }\mathbf{\varphi }=\mathbf{q}.$$

$${\psi }_m-{\psi }_n$$

$${\psi }_n$$

$${\varphi }_{\ell }$$

$$\mathrm{find}$$

$$s.to$$

$$\text{given}\{q_n{\}}_{n\in {\mathcal{N}}_q}\text{and}\{{\psi }_n{\}}_{n\in {\mathcal{N}}_{\psi }}.$$

$${\pi }_{\ell }^{mn}:=\{\begin{array}{ll}0& ,\text{if edge }\ell \notin {\mathcal{P}}_{mn}\\ +1& ,\text{if direction of }\ell \text{ agrees with path direction}\\ -1& ,\text{otherwise.}\end{array}$$

$$n_{\ell }^{\mathcal{C}}=\{\begin{array}{ll}0& ,\text{if edge }\ell \notin \mathcal{C}\\ +1& ,\text{if direction of }\ell \text{ agrees with cycle direction}\\ -1& ,\text{otherwise.}\end{array}$$

$$\mathrm{sign}({\mathbf{\varphi }}^{\prime }-\mathbf{\varphi })\odot {\mathbf{\pi }}^{mn}$$

$$\mathrm{sign}({\mathbf{\varphi }}^{\prime }-\mathbf{\varphi })\odot {\mathbf{\pi }}^{mn}>0$$

$$q_m^{\prime }>q_m\text{and}{q}_n^{\prime }<q_n$$

$$\underset{\mathbf{\varphi }}{\min }$$

$$s.to$$

$$\underset{\mathbf{\psi }}{\min }\frac{2}{3}\sum _{(m,n)\in \mathcal{P}}\frac{|{\psi }_m-{\psi }_n{|}^{\frac{3}{2}}}{\sqrt{a_{mn}}}-{\mathbf{q}}^{\top }\mathbf{\psi }.$$

$$\sum _{\ell =(m,n)\in \mathcal{P}}\mathrm{sign}({\mathbf{a}}_{\ell }^{\top }{\mathbf{\psi }}^{\ast })\sqrt{\frac{|{\mathbf{a}}_{\ell }^{\top }{\mathbf{\psi }}^{\ast }|}{a_{\ell }}}=q_n$$

$$(2x_{\ell }-1){\psi }_m+(1-2x_{\ell }){\psi }_n\ge a_{\ell }{\varphi }_{\ell }^2$$

$$x_{\ell }{\underline{\psi }}_n$$

$${\psi }_n+(x_{\ell }-1){\overline{\psi }}_n$$

$$2z_{\ell m}-2z_{\ell n}+{\psi }_n-{\psi }_m\ge a_{\ell }{\varphi }_{\ell }^2,$$

$$-{\overline{\varphi }}_{\ell }(1-x_{\ell })\le {\varphi }_{\ell }\le {\overline{\varphi }}_{\ell }{x}_{\ell },$$

$$\min$$

$$\mathrm{over}$$

$$s.to$$

$$r(\mathbf{\psi }):=\sum _{\begin{array}{c}(m,n)\in {\bar{\mathcal{P}}}_a\\ (m,n)\notin {\mathcal{S}}_{\mathcal{C}}^a\end{array}}|{\psi }_m-{\psi }_n|$$

$$\mathrm{sign}({\mathbf{\varphi }}_{\mathcal{C}}^{\prime }-{\mathbf{\varphi }}_{\mathcal{C}}+\lambda {\mathbf{n}}_{\mathcal{C}})\odot {\mathbf{n}}_{\mathcal{C}}=\mathrm{sign}\left((\lambda -{\lambda }_{\mathcal{C}}\right){\mathbf{n}}_{\mathcal{C}})\odot {\mathbf{n}}_{\mathcal{C}}>0.$$

$$G:=\underset{(m,n)\in {\bar{\mathcal{P}}}_a}{\max }\frac{|{\psi }_m-{\psi }_n|-a_{mn}{\varphi }_{mn}^2}{a_{mn}{\varphi }_{mn}^2}\ge 0$$

$${\mathbf{y}}_{t+1}={\mathbf{y}}_t-\mu [\mathbf{J}({\mathbf{y}}_t){]}^{-1}\mathbf{g}({\mathbf{y}}_t)$$

$$0$$

$$\stackrel{b}{=}{\pi }_{{\ell }_i}^{mn}\mathrm{sign}(w({\varphi }_{{\ell }_i}^{\prime })-w({\varphi }_{{\ell }_i}))$$

$$\stackrel{c}{=}\mathrm{sign}({\pi }_{{\ell }_i}^{mn})\mathrm{sign}(w({\varphi }_{{\ell }_i}^{\prime })-w({\varphi }_{{\ell }_i}))$$

$$\stackrel{d}{=}\mathrm{sign}({\pi }_{{\ell }_i}^{mn}w({\varphi }_{{\ell }_i}^{\prime })-{\pi }_{{\ell }_i}^{mn}w({\varphi }_{{\ell }_i}))$$

$$\stackrel{e}{=}\mathrm{sign}(({\psi }_{m_i}^{\prime }-{\psi }_{m_{i+1}}^{\prime })-({\psi }_{m_i}-{\psi }_{m_{i+1}})),$$

$${\psi }_{m_i}^{\prime }-{\psi }_{m_{i+1}}^{\prime }>{\psi }_{m_i}-{\psi }_{m_{i+1}}.$$

Full text

I Introduction

Natural gas has served as a critical energy source for decades, mainly for heating and electric power generation [1]. Thanks to the higher ramping capabilities of gas-fired generators, electric power system operators could achieve higher penetration of uncertain and intermittent renewable generation. In addition, the discovery of substantial new supplies of natural gas in the U.S. has led to a new thrust in development of gas-centered technologies and analytical tools [2].

Natural gas produced at gas pits and refineries is primarily transported to customer locations via a continent-wide network of pipelines [1]. The safe, reliable, and economical transportation of gas across these networks is ensured by gas system operators [3]. Considering the scale of natural gas networks (NGN), and their coupling with electric power grids, a plethora of analytical and computational challenges can be envisaged. Stand-alone and gas-electric coupled versions of network expansion planning, optimal scheduling, least-cost procurement, and security analysis are examples of problems that have gained increasing research interest; see [3], [4], [5], [6], [7]. These problems aim at optimizing varying objectives, while respecting network limitations and gas flow physics.

The flow of natural gas on pipelines is governed by partial differential equations, which under steady-state assumptions, yield nonlinear equations relating pressures and gas flows [8]. These equations reveal that the pressure drops along a pipe in the direction of flow due to friction. However, a minimum pressure needs to be maintained at consumer nodes to satisfy gas contracts. Therefore, compressors are placed on selected pipelines to increase the pressure at their output based on a typically multiplicative [1], and rarely additive law [9]. Operators need to solve the set of nonlinear equations governing gas flow in an NGN [10]: For each node, the operator fixes the gas pressure or gas injection rate to specified values. Given also the compression ratios, the GF task aims at finding the injections and pressures at all nodes, as well as the gas flows on all pipes. While solving the GF task is central for numerous NGN operations, it is hard to do so even under steady-state and balanced conditions for non-tree networks [1].

The GF task is usually handled by Newton-Raphson (NR)-based solvers. However, their convergence can be sensitive to initialization [11]. A semidefinite program (SDP)-based GF solver attaining a higher success probability than the NR scheme, is developed in [10]. Nevertheless, the SDP based solver fails to solve the GF problem if the network state is far from the states considered in designing the solver. The necessity of proper initialization may be avoided for simpler networks without compressors as the flows and pressures may be found as optimal primal-dual solutions of a convex minimization [12]. Nevertheless, for practical meshed NGNs with compressors, an initialization-independent GF solver is still a research pursuit [8]. Setting scalability aside, if one uses a nonlinear solver for the GF task, the uniqueness of a solution becomes critical. References [9] and [13] prove the uniqueness of a GF solution for NGNs with additive compressors.

The contribution of this work is on four fronts: First, Section III establishes that the nonlinear steady-state GF equations enjoy a unique solution even with multiplicative compressors. Building on [8] where uniqueness was shown for GF setups with a single fixed-pressure node, here uniqueness is non-trivially generalized to setups with multiple fixed-pressure nodes. Second, Section IV reformulates the GF task as a convex minimization. The obtained solver can handle GF setups with a single fixed-pressure node and compressors not on cycles. Third, Section V expands the analytical claims for the MI-QCQP gas flow solver of [8]. Different from the convex minimization approach, this solver applies to any GF setup and any network. The MI-QCQP solver introduces binary variables to capture flow directions; relaxes the nonlinear GF equations to quadratic inequalities; and uses a carefully selected objective to promote the exactness of the relaxation. The relaxation is provably exact in NGNs with non-overlapping cycles and a single fixed-pressure node. This significantly extends the claim of [8], where exactness was proved for non-overlapping cycles and a single fixed-pressure node, but did not allow for compressors in cycles. Having compressors in cycles is a typical arrangement, e.g., when two compressors are connected in parallel. Fourth, to accelerate the MI-QCQP solver, the bilinear terms involved are handled through McCormick linearization. Numerical tests on meshed networks with overlapping cycles and multiple fixed-pressure nodes demonstrate that the MI-QCQP solver finds the unique GF solution even when the assumed sufficient conditions are violated.

II Gas Flow Problem

A natural gas network (NGN) can be represented by a directed graph $\mathcal{G}=(\mathcal{N},\mathcal{P})$. The nodes in the graph represent points of gas supply, demand, or network junctions. The edges are directed, and represent pipelines or compressors. Nodes are indexed by $n\in\mathcal{N}:=\{1,\cdots,N\}$ and edges by $\ell\in\mathcal{P}:=\{1,\ldots,P\}$. Each edge $\ell=(m,n)$ is assigned a direction from the origin node $m$ to the destination node $n$. If $(m,n)\in\mathcal{P}$, then $(n,m)\notin\mathcal{P}$. For edges corresponding to pipes, this direction is selected arbitrarily. For edges denoting compressors, the direction coincides with the direction of gas flow, since compressors allow only unidirectional flow of gas.

For each node $n\in\mathcal{N}$, let $q_{n}$ be the gas injection rate from node $n$ to the NGN. By convention, the gas injection $q_{n}$ is positive for gas source nodes; negative for demand nodes; and zero for junction nodes. Vector $\mathbf{q}\in\mathbb{R}^{N}$ collects the gas injections across all nodes.

For each edge $\ell=(m,n)\in\mathcal{P}$, let $\phi_{\ell}$ denote its gas flow rate. By convention, the flow $\phi_{\ell}$ is positive if gas flows from node $m$ to $n$; and negative, otherwise. The conservation of mass at each node $n\in\mathcal{N}$ dictates that

[TABLE]

Under steady-state conditions, the input and output flows on a pipe are identical, and so gas injections are balanced at all times, that is $\sum_{n=1}^{N}q_{n}=0$. Because of this, from the $N$ linear equations in (1), only $(N-1)$ are linearly independent.

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

[TABLE]

Using $\mathbf{A}$, equation (1) can be compactly expressed as

[TABLE]

where vector $\boldsymbol{\phi}\in\mathbb{R}^{P}$ stacks the flows $\phi_{\ell}$’s along all edges.

For medium- and high-pressure networks, the gas flows on pipelines relate to nodal pressures through a set of nonlinear partial differential equations [14], [15]. These equations model the gas flow dynamics evolving across time and spatially along the pipeline length. However, simplifying assumptions such as ignoring friction, geographical tilt, variations in ambient temperature, and time-varying gas injections, yield the popular steady-state Weymouth equation [16]. If $\psi_{n}>0$ denotes the squared gas pressure at node $n\in\mathcal{N}$, the pressure drop across pipeline $\ell=(m,n)\in\mathcal{P}$ is given by

[TABLE]

where parameter $a_{\ell}>0$ depends on physical properties of the pipeline [14]. The function $\operatorname{sign}(x)$ returns $+1$ if $x>0$; $-1$ if $x<0$; and [math] if $x=0$. The absolute value in (3a) signifies that pressure drops along the direction of flow. In particular, the drop in squared pressures is proportional to the squared flow. We will henceforth refer to $\psi_{m}$ as pressure rather than squared pressure for brevity. Let us collect all $\psi_{n}$’s in $\boldsymbol{\psi}\in\mathbb{R}^{N}$.

To enable the desired flow of gas in an NGN while maintaining pressures within acceptable limits, system operators install compressors at selected pipelines. A pipeline hosting a compressor can be modeled by an ideal compressor which increases the gas pressure, followed by a lossy pipeline that incurs a pressure drop per (3). Apparently, the gas flows on the two edges are identical. Let the subset of edges hosting ideal compressors be $\mathcal{P}_{a}\subset\mathcal{P}$. The edges in $\mathcal{P}_{a}$ are also referred to as active pipelines. The pressures across an active pipeline or compressor $\ell=(m,n)\in\mathcal{P}_{a}$ are related as

[TABLE]

where $\alpha_{\ell}>0$ is the multiplicative compression factor for compressor $\ell$. The unidirectional flow permitted for a compressor is enforced by (4b). The remaining edges, that is the edges not hosting ideal compressors, constitute the set $\bar{\mathcal{P}}_{a}:=\mathcal{P}\setminus\mathcal{P}_{a}$ and abide by (3) instead of (4).

In an NGN, a node $r\in\mathcal{N}$ is selected as a reference node. Its pressure is kept fixed. Given $\psi_{r}$, if the nodal pressures $\boldsymbol{\psi}$ are known, the flows $\boldsymbol{\phi}$ can be readily computed; and vice versa. This fact follows immediately from (3)–(4), and is itemized as the next lemma to be used in subsequent arguments.

Lemma 1**.**

Given a reference pressure $\psi_{r}$ for some $r\in\mathcal{N}$, a pair $(\boldsymbol{\phi},\boldsymbol{\psi})$ satisfying (3)–(4) is uniquely characterized by either $\boldsymbol{\phi}$ or $\boldsymbol{\psi}$.

The task of finding $\boldsymbol{\phi}$ or $\boldsymbol{\psi}$ given a combination of nodal injections and pressures constitutes the gas flow (GF) problem. Oftentimes, gas supply nodes are tuned to maintain a fixed pressure while injecting variable amounts of gas to meet the prescribed pressure under variable demands [3], [17]. Let set $\mathcal{N}_{\psi}\subset\mathcal{N}$ consist of all nodes with fixed pressures $\psi_{n}$’s. The reference node $r$ belongs to $\mathcal{N}_{\psi}$ by definition. Its complement set $\mathcal{N}_{q}:=\mathcal{N}\setminus\mathcal{N}_{\psi}$ consists of all nodes with fixed injections $q_{n}$’s. Then, the GF problem can be formally stated now.

Definition 1**.**

Given pressures $\psi_{n}$ for $n\in\mathcal{N}_{\psi}$; injections $q_{n}$ for $n\in\mathcal{N}_{q}$; the ratios $\alpha_{\ell}$ for all compressors $\ell\in\mathcal{P}_{a}$; and the friction parameters $a_{\ell}$ for all lossy pipes $\ell\in\bar{\mathcal{P}}_{a}$, the GF problem aims at finding the triplet $(\boldsymbol{\psi},\boldsymbol{\phi},\mathbf{q})$ satisfying the GF equations (2)–(4).

The GF task involves $N-1+P$ equations over $N-1+P$ unknowns. It can be posed as the feasibility problem

[TABLE]

Albeit (2) and (4) are linear, the piecewise quadratic Weymouth equation in (3) is non-convex, while the requirement $\{\phi_{\ell}\geq 0\}_{\ell\in\mathcal{P}_{a}}$ further complicates the task. The GF problem is typically solved using the Newton-Raphson’s method, yet its convergence depends on the initialization [11], [17], [10]. Commercially available software require careful manual tuning by gas network operator personnel, though that could be attributed to more detailed models of NGN components.

A popular rendition of the GF problem considers the reference node as the only fixed-pressure node, and all other nodes as fixed-injection nodes [18], [10], [8]. For this rendition, solving the GF problem becomes trivial for a tree network by inverting (2) and using Lemma 1. However, for a meshed NGN, solving the GF problem remains non-trivial. Before developing new GF solvers, the next section establishes that the GF problem in (G1) enjoys a unique solution.

III Uniqueness of the GF Solution

We commence with the uniqueness of the GF task under the setup of a single fixed-pressure node, proved in [8, Th. 1].

Theorem 1** ([8]).**

If $\mathcal{N}_{\psi}=\{r\}$ and $\mathcal{N}_{q}=\mathcal{N}\setminus\{r\}$, the gas flow problem (G1) has a unique solution, if feasible.

Although the single fixed-pressure setup has been studied widely, setups with multiple fixed-pressure nodes are of critical interest too. This is because gas is typically injected at supplier sites using a controller that maintains constant pressure, rather than constant rate. To address this need, this section builds on Th. 1 and establishes the uniqueness of the steady-state GF equations for any $(\mathcal{N}_{\psi},\mathcal{N}_{q})$ setup. Before doing so, let us briefly review some graph theory preliminaries.

A directed graph $\mathcal{G}=(\mathcal{N},\mathcal{P})$ is connected if there exists a sequence of adjacent edges between any two nodes. All graphs considered in this work are assumed to be connected. A sequence of adjacent edges between nodes $m$ and $n$ constitutes a path $\mathcal{P}_{mn}\subset\mathcal{P}$. The directionality assigned to path $\mathcal{P}_{mn}$ is from $m$ to $n$. Note that nodes $m$ and $n$ could be connected by multiple paths. Thus, with slight abuse in notation, path $\mathcal{P}_{mn}$ shall represent any arbitrary path between $m$ and $n$, unless additional conditions are provided. For path $\mathcal{P}_{mn}$, we can define an indicator vector $\boldsymbol{\pi}^{mn}\in\{0,1\}^{P}$ with $\ell-$th entry

[TABLE]

A cycle is a sequence of adjacent edges (without edge or node repetition) that starts and ends at the same node. With a slight abuse of terminology, the statement ‘cycle $\mathcal{C}$ contains node $i$’ will mean that there exists an edge in $\mathcal{C}$ that is incident to node $i$. For any cycle $\mathcal{C}$, we can select an arbitrary direction and define its indicator vector $\mathbf{n}^{\mathcal{C}}$ with $\ell$-th entry

[TABLE]

A tree is a connected graph with no cycles.

After the graph theoretic preliminaries, we proceed with the uniqueness of the GF solution for the general GF setup. This proof builds upon the ensuing two lemmas, which are proved in the appendix.

Lemma 2**.**

Consider path $\mathcal{P}_{mn}$ along edges $\{\ell_{1},\dots,\ell_{k}\}$ with indicator $\boldsymbol{\pi}^{mn}$. For fixed pressures $\psi_{m}$ and $\psi_{n}$, if flow vectors $\boldsymbol{\phi}$ and $\boldsymbol{\phi}^{\prime}$ with $\boldsymbol{\phi}\neq\boldsymbol{\phi}^{\prime}$, satisfy (3)–(4), they cannot satisfy

[TABLE]

where the strict inequalities are understood entrywise.

To get some intuition, suppose that $\boldsymbol{\pi}^{mn}$ takes the value of $+1$ for edges $\{\ell_{1},\dots,\ell_{k}\}$, and [math] for the remaining edges. According to Lemma 2, if two pairs $(\boldsymbol{\phi},\boldsymbol{\psi})$ and $(\boldsymbol{\phi}^{\prime},\boldsymbol{\psi}^{\prime})$ satisfy (3)–(4) with $\psi_{m}=\psi_{m}^{\prime}$ and $\psi_{n}=\psi_{n}^{\prime}$, then the flows along $\mathcal{P}_{mn}$ cannot uniformly increase from $\boldsymbol{\phi}$ to $\boldsymbol{\phi}^{\prime}$. In other words, $\phi_{\ell}^{\prime}>\phi_{\ell}$ cannot occur simultaneously for all $\ell\in\mathcal{P}_{mn}$. Flows cannot uniformly decrease either ($\phi_{\ell}^{\prime}<\phi_{\ell}$ for all $\ell\in\mathcal{P}_{mn}$). This holds merely because the pressure drop across a pipe decreases monotonically with gas flow from and compressors perform a linear scaling [cf. (3) and (4)].

The next lemma describes an interesting effect on how gas flows get redistributed when gas injections change.

Lemma 3**.**

Consider two pairs $(\mathbf{q},\boldsymbol{\phi})$ and $(\mathbf{q}^{\prime},\boldsymbol{\phi}^{\prime})$ satisfying (2). If $\mathbf{q}\neq\mathbf{q}^{\prime}$, there exists a path $\mathcal{P}_{mn}$ between nodes $m$ and $n$ such that

[TABLE]

where $\boldsymbol{\pi}^{mn}$ is the indicator vector for $\mathcal{P}_{mn}$.

Lemma 3 predicates that if gas injections change, there exists a path: i) along which flows increase uniformly; ii) the source node of the path has increased injection; and iii) the destination node has decreased injection. Lemma 3 has been established in [9] via mathematical induction; see the appendix for an alternative perhaps more intuitive proof.

Using Theorem 1 and Lemmas 2–3, we next prove the uniqueness of the GF task under the general setup.

Theorem 2**.**

The gas flow problem (G1) has a unique solution, if feasible.

Proof.

Proving by contradiction, assume $(\mathbf{q},\boldsymbol{\phi},\boldsymbol{\psi})$ and $(\mathbf{q}^{\prime},\boldsymbol{\phi}^{\prime},\boldsymbol{\psi}^{\prime})$ are two distinct solutions of (G1). Consider the GF setup where $|\mathcal{N}_{\psi}|>1$; the special case of $|\mathcal{N}_{\psi}|=1$ is covered by Theorem 1. If $\mathbf{q}\neq\mathbf{q}^{\prime}$, then Lemma 3 implies that there exists a path $\mathcal{P}_{mn}$ with indicator vector $\boldsymbol{\pi}^{mn}$ satisfying (6a). Moreover, it holds that $q_{m}^{\prime}>q_{m}$ and $q_{n}^{\prime}<q_{n}$ from (6b). By definition, gas injections $q_{i}$ are fixed for all nodes $i\in\mathcal{N}_{q}$. Therefore, nodes $m$ and $n$ cannot be fixed-injection nodes. They have to be fixed-pressure nodes belonging to $\mathcal{N}_{\psi}$, implying $\psi_{m}=\psi_{m}^{\prime}$ and $\psi_{n}=\psi_{n}$. However, with the pressures at nodes $m$ and $n$ fixed, the inequality (6a) contradicts Lemma 2. Hence, the assumption of unequal injections is refuted, implying $\mathbf{q}=\mathbf{q}^{\prime}$.

Given $\mathbf{q}$ and the reference pressure $\psi_{r}$, Theorem 1 asserts that there is unique triplet $(\mathbf{q},\boldsymbol{\phi},\boldsymbol{\psi})$ satisfying the GF equations. Since $\mathbf{q}=\mathbf{q}^{\prime}$, the triplets $(\mathbf{q},\boldsymbol{\phi},\boldsymbol{\psi})$ and $(\mathbf{q}^{\prime},\boldsymbol{\phi}^{\prime},\boldsymbol{\psi}^{\prime})$ have to coincide, which completes the proof. ∎

The uniqueness claim of Theorem 2 is fairly general, since it applies to any NGN topology and any GF setup with a single or multiple fixed-pressure nodes. Having established uniqueness, the next two sections develop a suite of GF solvers: Section IV builds upon an existing convex solver for GF setups with a single fixed-pressure node and no compressors. We develop an unconstrained convex solver as well as an extension that handles compressors on non-overlapping cycles. Section V adopts a convex relaxation and puts forth an MI-QCQP to handle more general GF setups. The relaxation is provably exact for NGNs with a single fixed-pressure node and non-overlapping cycles. Nonetheless, numerical tests demonstrate that this MI-QCQP succeeds in finding the unique GF solution in NGNs with multiple fixed-pressure nodes and overlapping cycles as long as compressors are not on overlapping cycles.

IV Energy Function Minimization

This section studies the GF task for the special case of $|\mathcal{N}_{\psi}|=1$. In an NGN without compressors, the GF task is posed as a convex minimization. The approach can be extended to networks having compressors, but not on cycles.

IV-A Existing Constrained Energy Function-based GF Solver

Consider solving the GF task for a single fixed-pressure node (the reference node $r$) and in an NGN without compressors. This task boils down to solving equations (2)–(3). As shown in [12], the gas flows $\boldsymbol{\phi}$ for this GF setup can be found as the minimizer of the convex minimization

[TABLE]

This can be readily verified by the first-order optimality conditions of (7). In addition, the pressures $\boldsymbol{\psi}$ can be recovered from the optimal Lagrange multipliers $\boldsymbol{\xi}\in\mathbb{R}^{N}$ associated with constraint (7b): If $\boldsymbol{\xi}$ is shifted by a constant so that its $r$-th entry equals $\psi_{r}$, the remaining entries of this shifted $\boldsymbol{\xi}$ equal $\boldsymbol{\psi}$. Problem (7) can be reformulated as a second-order cone program or tackled via dual decomposition; see [19].

IV-B Novel Unconstrained Energy Function-based GF Solver

Rather than solving (7) over $\boldsymbol{\phi}$, here we show that one can alternatively find the GF solution via an unconstrained convex minimization over $\boldsymbol{\psi}$ as

[TABLE]

The convexity of this objective function follows from composition rules. Since this function is convex and differentiable, its unconstrained minimization is equivalent to nulling its gradient vector. Setting the $n$-th entry of this gradient to zero reveals that the minimizer $\boldsymbol{\psi}^{*}$ of (8) satisfies

[TABLE]

where $\mathbf{a}_{\ell}^{\top}$ is the $\ell$-th row of matrix $\mathbf{A}$. Equation (9) is equivalent to eliminating the flows $\boldsymbol{\phi}$ from (2) and (3). As with (7), the ambiguity in pressures could be handled by shifting $\boldsymbol{\psi}^{*}$ by a constant, so that $\psi_{r}^{*}$ agrees with the given pressure at the reference node $r$. Once pressures $\boldsymbol{\psi}^{*}$ have been determined, flows can be found using Lemma 1.

Remark 1**.**

In the absence of compressors and when $|\mathcal{N}_{\psi}|=1$, the GF task becomes structurally similar to the water flow problem in water distribution networks without pumps [19]. Therefore, the (un)-constrained energy function minimization approaches of (7)–(8) apply to the gas flow and water flow problems alike. For water networks, the decomposition technique of [19] extends (7)–(8) to water network setups with $|\mathcal{N}_{\psi}|=1$ and pumps, but pumps cannot lie on cycles. A similar technique can be used to solve the GF problem with compressors not on cycles and $|\mathcal{N}_{\psi}|=1$. The only modification needed relates to accounting for the multiplicative pressure law in gas compressors [cf. (4a)] vis-à-vis the additive pressure law of water pumps. Additionally, the decomposition algorithm may be extended to accommodate compressors on non-overlapping cycles using the flow-recovery procedure provided later as Algorithm 1. Since carrying over this decomposition technique from the water flow to the gas flow context is straight-forward and due to space limitations, it is not presented here.

To handle GF setups with $|\mathcal{N}_{\psi}|>1$ and/or NGNs with compressors in loops, a convex relaxation of the Weymouth equation is pursued in the next section.

V MI-QCQP Relaxation

The minimization approaches of (7)–(8) provide computationally efficient methods to solve the GF problem, but exhibit three limitations: i) they cannot handle multiple fixed-pressure nodes ($|\mathcal{N}_{\psi}|>1$); ii) cannot handle compressors on cycles; and iii) cannot be extended to optimal gas flow formulations (e.g., along the lines of [20]). To overcome these limitations, this section presents an MI-QCQP-based solver that is applicable to any GF setup.

V-A Problem Reformulation

The non-convexity of (G1) is due to the Weymouth equation in (3a). The piecewise quadratic equalities can be relaxed to convex inequality constraints: The pressure drop along a lossy pipe $\ell=(m,n)\in\bar{\mathcal{P}}_{a}$ is relaxed to

• •

$\psi_{m}-\psi_{n}\geq a_{\ell}\phi_{\ell}^{2}$ for $\phi_{\ell}\geq 0$; or

• •

$\psi_{n}-\psi_{m}\geq a_{\ell}\phi_{\ell}^{2}$ for $\phi_{\ell}\leq 0$.

The two cases can be differentiated using a binary variable $x_{\ell}$ capturing the direction of flow $\phi_{\ell}$. The relaxed pressure drop equations can be compactly written as

[TABLE]

where $x_{\ell}=1$ corresponds to $\phi_{\ell}\geq 0$; and $x_{\ell}=0$ to $\phi_{\ell}\leq 0$. Despite the relaxation, the bilinear terms $x_{\ell}\psi_{m}$ make the aforementioned constraint non-convex.

The McCormick linearization, popular for approximating multilinear terms by their linear convex envelopes, can be used to handle these bilinear terms [21]. For the special case of bilinear terms involving at least one binary term, the McCormick linearization becomes exact. In fact, it is related to the so termed big-$M$ trick, but instead of using a single arbitrarily large value for $M$, it selects different values for $M$ that are specialized per product of variables, which could potentially reduce the running time of mixed-integer programming solvers. Let us briefly review the linearization. Consider the constraint $z_{\ell n}=x_{\ell}\psi_{n}$, for which $x_{\ell}\in\{0,1\}$ and $\psi_{n}\in[\underline{\psi}_{m},\overline{\psi}_{n}]$. This constraint can be equivalently expressed via four linear inequalities

[TABLE]

To verify the exactness, observe that when $x_{\ell}=1$, constraint (10b) yields $z_{\ell n}=\psi_{n}$ and (10a) holds trivially. When $x_{\ell}=0$, constraint (10a) enforces $z_{\ell n}=0$ and (10b) holds trivially. Hence, the constraints in (10) ensure that $z_{\ell n}=x_{\ell}\psi_{n}$.

To arrive at an MI-QCQP relaxation of (G1), for all lossy pipes $\ell\in\bar{\mathcal{P}}_{a}$, the pressure drop constraint of (3a) is replaced by (10) and

[TABLE]

where $\overline{\phi}_{\ell}$ is an upper bound on $|\phi_{\ell}|$. Constraint (11a) represents the relaxed Weymouth equation, and constraint (11b) defines $x_{\ell}=\operatorname{sign}(\phi_{\ell})$. Similar relaxations have been previously used in [4], [22], [8]; see Section VI for a detailed comparison.

When solving the GF problem with the Weymouth equations relaxed, the obtained solution is useful only if the relaxation is exact, that is when (11a) holds with equality for all $\ell$. To render the relaxation provably exact, we convert the feasibility problem (G1) to the MI-QCQP minimization

[TABLE]

The optimization variable $\mathbf{x}$ stacks $\{x_{\ell}\}_{\ell\in\bar{\mathcal{P}}_{a}}$, and the objective function is judiciously selected as

[TABLE]

where $\mathcal{S}_{\mathcal{C}}^{a}$ is the set of cycles with compressors. These cycles will be also termed as active cycles. The cost $r(\boldsymbol{\psi})$ sums up the absolute pressure differences across all lossy pipes not in active cycles. Despite the non-convexity of (G2) due to the binary variables, this minimization can be handled for moderately sized networks thanks to the advancements in mixed-integer second-order cone solvers. The computational performance of (G2) is further corroborated by our tests. The next section provides network conditions under which the exactness of (G2) can be guaranteed analytically. The tests in Section VII demonstrate numerically that solving (G2) renders the relaxation exact for a much broader class of networks.

For solving tasks such as (G1), NR-based or fixed-point iteration solvers are often preferred as opposed to optimization-based solvers due to computational superiority. However, in addition to guaranteeing convergence irrespective of initialization, problem (G2) can also be used as follows:

• •

Infeasibility: As a relaxation of (G1), (G2) can be used to screen infeasible GF instances; see Section VII for tests. Such screening is of practical use as suggested in [23].

• •

Initialization: Problem (G2) could be terminated before reaching optimality to yield initializations for NR solvers, hence combining the benefits of both approaches.

• •

Optimal gas flow: The cost of (G2) could be useful as a penalty term that can be added to optimization problems [19], [20]. However, guaranteeing exact relaxation for such problems would need further analysis.

V-B Exactness of the Relaxation

The relaxation in (G2) will be analytically shown to be exact under the following network conditions.

Condition 1**.**

The GF setup has a single fixed-pressure node, that is $|\mathcal{N}_{\psi}|=1$.

Condition 2**.**

Each edge of the NGN belongs to at most one cycle.

Condition 3**.**

The NGN does not exhibit circulation of gas, that is $\mathbf{n}^{c}\odot\boldsymbol{\phi}\not>\mathbf{0}$ and $\mathbf{n}^{c}\odot\boldsymbol{\phi}\not<\mathbf{0}$ for every cycle $\mathcal{C}$.

Under Condition 1, the nodal injections are fixed a priori and the GF task aims at finding the associated $(\boldsymbol{\psi},\boldsymbol{\phi})$. Albeit Definition 1 considered the GF task with multiple fixed-pressure nodes, the setup of a single fixed-pressure node is commonly met; see [18], [10], [8], [24]. Regarding Condition 2, although it may seem restrictive at the outset, it is satisfied by several practical gas networks [25]. For Condition 3, a circulation occurs when gas flows around a cycle along the same direction. It is easy to verify that gas cannot circulate in a cycle without compressors, since the incurred pressure drops along the cycle will all be in the same direction and thus cannot sum up to zero. In cycles with compressors, gas circulation can occur though it would cause an undesirable loss of energy. However, the tests of Section VII demonstrate that the relaxation in (G2) is exact even in setups where the sufficient Conditions 1–3 are all violated.

The next exactness claim applies to the GF setup with known injections. From Lemma 1, we know that solving the GF task is equivalent to finding the correct flows $\boldsymbol{\phi}$. The next result provides conditions under which (G2) yields flows $\boldsymbol{\phi}$ with partially correct entries. An algorithm to retrieve the entire $\boldsymbol{\phi}$ and thus eventually solve (G1) is presented afterwards.

Theorem 3**.**

Let $\boldsymbol{\phi}$ be the unique flow vector solving (G1), and $\boldsymbol{\phi}^{\prime}$ the flow vector minimizing (G2). Under Cond. 1–3, it holds $\phi_{\ell}^{\prime}=\phi_{\ell}$ for all edges $\ell$ not belonging to active cycles.

Theorem 3 establishes that the only possible mismatches between $\boldsymbol{\phi}^{\prime}$ and $\boldsymbol{\phi}$ occur only at the edges lying on cycles with compressors. Then, if there are no cycles with compressors, the GF problem is solved correctly; see also [8, Th. 2].

Corollary 1**.**

Under Cond. 1–2, for an NGN without compressors in cycles, the minimizer of (G2) solves (G1) as well.

Corollary 1 identifies a setup where (G2) is equivalent to solving (G1). Nonetheless, if there are no compressors in cycles, one would prefer tackling (G1) using the solvers of Section IV. This is because running the decomposition technique discussed in Remark 1 and solving (7), are simpler than solving the MI-QCQP of (G2).

V-C Recovering the GF Solution

Returning to the general setup, we next provide a procedure to retrieve the solution $\boldsymbol{\phi}$ of (G1) given a minimizer $\boldsymbol{\phi}^{\prime}$ of (G2). From Theorem 3, vector $\boldsymbol{\phi}^{\prime}$ needs to be corrected only at the entries corresponding to edges in active cycles. To this end, we first put forth an algorithm to correct the flows within a single active cycle, and then delineate the steps to systematically correct the flows for all active cycles of the network.

Consider an active cycle $\mathcal{C}$ with $N_{\mathcal{C}}$ nodes. Let $\psi_{0}$ be a known pressure on node $0\in\mathcal{C}$, and $\boldsymbol{\phi}_{\mathcal{C}}^{\prime}$ be the $N_{\mathcal{C}}$-length subvector of $\boldsymbol{\phi}^{\prime}$ collecting the flows on $\mathcal{C}$. Similarly, let $\mathbf{n}_{\mathcal{C}}$ be the $N_{\mathcal{C}}$-length subvector of the indicator vector for cycle $\mathcal{C}$. The next lemma explains how $\boldsymbol{\phi}_{\mathcal{C}}$ can be recovered from $\boldsymbol{\phi}_{\mathcal{C}}^{\prime}$.

Lemma 4**.**

Given a known pressure $\psi_{0}$, and flows $\boldsymbol{\phi}_{\mathcal{C}}^{\prime}$ on active cycle $\mathcal{C}$ obtained from (G2), Algorithm 1 determines the corrected gas flows $\boldsymbol{\phi}_{\mathcal{C}}$ such that the relaxed Weymouth equations in (11) are satisfied with equality.

Proof.

Because $\boldsymbol{\phi}$ and $\boldsymbol{\phi}^{\prime}$ both satisfy (2), it follows that $(\boldsymbol{\phi}-\boldsymbol{\phi}^{\prime})\in\operatorname{null}(\mathbf{A}^{\top})$. Since there are no overlapping cycles, we have that $\boldsymbol{\phi}_{\mathcal{C}}=\boldsymbol{\phi}_{\mathcal{C}}^{\prime}+\lambda_{\mathcal{C}}\mathbf{n}_{\mathcal{C}}$ for some $\lambda_{\mathcal{C}}\in\mathbb{R}$. To recover $\boldsymbol{\phi}_{\mathcal{C}}$, we next provide a method for finding $\lambda_{\mathcal{C}}$.

Suppose one is given a $\lambda\in\mathbb{R}$. Given pressure $\psi_{0}$ and the candidate flow vector $\boldsymbol{\phi}_{\mathcal{C}}^{\prime}+\lambda\mathbf{n}_{\mathcal{C}}$, one can calculate the pressures along $\mathcal{C}$ sequentially using (3a) and (4a). Upon completing the cycle, the pressure at node $0\in\mathcal{C}$ will be evaluated to the value of $\hat{\psi}_{0}(\lambda)$. The value $\hat{\psi}_{0}(\lambda)$ may not be equal to $\psi_{0}$. Note that for $\lambda>\lambda_{\mathcal{C}}$, it holds that

[TABLE]

Using the above along with the argument used in the proof of Lemma 2, it can be shown that $\hat{\psi}_{0}(\lambda)<\psi_{0}$. In a similar fashion, if $\lambda<\lambda_{\mathcal{C}}$, then $\hat{\psi}_{0}(\lambda)>\psi_{0}$. Therefore, the function $\hat{\psi}_{0}(\lambda)-\psi_{0}$ is monotonic in $\lambda$, and $\hat{\psi}_{0}(\lambda)=\psi_{0}$ if and only if $\lambda=\lambda_{\mathcal{C}}$. Thanks to this monotonicity, one can find $\lambda_{\mathcal{C}}$ iteratively using bisection, tabulated as Algorithm 1. ∎

Lemma 4 shows that $\boldsymbol{\phi}_{\mathcal{C}}$ can be recovered from $\psi_{0}$ and $\boldsymbol{\phi}_{\mathcal{C}}^{\prime}$ using a bisection technique on $\lambda$. The limits for the search space $[\underline{\lambda},\overline{\lambda}]$ of $\lambda$ can be found using engineering constraints on gas flows. In fact, these limits can be tightened since the entries of $\boldsymbol{\phi}_{\mathcal{C}}^{\prime}$ and $\boldsymbol{\phi}_{\mathcal{C}}(\lambda)=\boldsymbol{\phi}_{\mathcal{C}}^{\prime}+\lambda\mathbf{n}_{\mathcal{C}}$ corresponding to any compressor in $\mathcal{C}$ must have the same sign due to (4b).

We next provide the steps to find the correct GF solution using the flow $\boldsymbol{\phi}^{\prime}$ obtained from (G2):

T1) Select a spanning tree $\mathcal{T}$ of the NGN graph $\mathcal{G}$ rooted at the reference node $r$.

T2) Starting from node $r$, traverse $\mathcal{T}$ via a depth-first search.

T3) If a node $n$ does not belong to an active cycle of $\mathcal{G}$, calculate its pressure as follows: If the edge connecting node $n$ to its parent node in $\mathcal{T}$ is a lossy pipe, use (3a); else, if this edge is a compressor, use (4a).

T3) If a node $n$ belongs to an active cycle $\mathcal{C}$, check if the flows in cycle $\mathcal{C}$ have been corrected. If the flows are already corrected or if $i$ is the first node in $\mathcal{C}$ that is encountered, compute the nodal pressure as in step T3). Else, pass the pressure at the parent node of $i$ (which is also in $\mathcal{C}$) along with the non-corrected flow subvector $\boldsymbol{\phi}_{\mathcal{C}}^{\prime}$ to Algorithm 1 and obtain the corrected flows on $\mathcal{C}$.

T4) Continue until all nodes in $\mathcal{T}$ have been traversed.

VI Comparison to Prior Work

This work puts forth three novel components: c1) proving the uniqueness of the GF problem solution under steady-state conditions; c2) proposing GF solvers based on energy function minimization; and c3) devising a provably exact MI-QCQP relaxation. These components are next contrasted to existing related works:

c1) Uniqueness: For an NGN with no compressors, the GF solution may be found as a minimizer of (7); see [12], [23]. A linearization technique has also been put forth to accelerate solving (7) [23]. References [9] and [13] broaden the uniqueness claim for NGNs with additive compressors of constant gain. These works formulate strictly convex problems that yield a GF solution; hence proving uniqueness by convexity. However, gas compressors are oftentimes multiplicative, so that the previous uniqueness claims do not carry over. In our work [8], uniqueness was proved for multiplicative compressors under any network topology, but with a single fixed-pressure node. Theorem 2 generalizes all past claims for NGNs with multiplicative compressors, any topology, and an arbitrary number of fixed-pressure nodes.

c2) Energy Function Minimization: Problem (7) dates back to [12], and has since been used for solving the GF task; verifying the feasibility and uniqueness of a GF instance [23]; and initializing optimization problems. However, its applicability was limited to NGNs without compressors. As explained in Remark 1, this work suggests using (7) to handle NGNs with compressors on non-overlapping cycles. We also present the unconstrained energy function formulation of (8).

c3) MI-QCQP relaxation of GF: The key difficulty in solving (optimal) GF problems stems from the non-linear Weymouth equation. A disjunctive convex relaxation of this equation was found to be efficient in [4], [22]. Numerous studies have thereon employed similar convex relaxations; see [5], [7], [26]. Unfortunately, it is hard to guarantee the exactness of these relaxations. An effective heuristic is to fix the binary variables involved to the values obtained by the convex relaxation and handle the resultant non-convex nonlinear program through a general solver [4], [5]. A gas-electric flow problem was solved in [26], wherein a cost function was proposed that was numerically found useful towards attaining exact relaxation. Unlike previous works, the MI-QCQP formulation of Section V provides theoretical guarantees for exact relaxation, while expanding the claims of [8]. The GF solver developed in [8] was applicable to NGN’s with compressors not on cycles. However, in this work, the cost function of (G2) is meticulously designed to ensure that correct flows are obtained outside active cycles. Additionally, Algorithm 1 is developed to enable flow correction on active cycles efficiently. Although the GF problem is intrinsically simpler than the optimal gas (and possible electric) flow problem considered in prior works, this work lays a foundation towards analytical guarantees for exact relaxation. It has been recently shown that exact relaxation of network flow optimization problems may be guaranteed using a convex penalty [27]. It is worth mentioning that a related MI-QCQP formulation of the water flow problem in [19], can also provably yield an exact convex relaxation for the optimal water flow task [20].

VII Numerical Tests

The proposed GF solver based on the relaxed MI-QCQP (G2) and Algorithm 1 was tested on the modified Belgian benchmark NGN and the GasLib-40 NGN of Fig. 1. Starting with the Belgian NGN, the pipe coefficients and compressor ratios were derived based on the nodal pressures and edge flows reported in [12]. The network contains three compressors, which are modeled as ideal compressors followed by lossy pipes. Problem (G2) was solved using the MATLAB-based optimization toolbox YALMIP using CPLEX as the MI-QCQP solver [28], [29]. All tests were conducted on a 2.7 GHz Intel Core i5 computer with 8 GB RAM.

As a model validation step, we first tested the (G2) solver on the original Belgian network, which is a tree, except for one cycle formed by parallel compressors, see Fig. 1. The pressure at node $1$ was treated as reference. The flow values obtained from (G2) agreed with those of [12] for all edges except for the edges along the active cycle. Similarly, the pressures agreed for all nodes other than node $20$. Therefore, the pressure at node $19$ and the flows on edges $(19,21)$, $(19,22)$, $(20,21)$, $(20,22)$ were passed to Algorithm 1 for correction. The final result was found to coincide with [12].

The Belgian network was subsequently augmented by additional pipelines; see Fig. 1. The resulting modified network has overlapping cycles, thus violating Condition 2 required in Theorem 3. To get reasonable friction coefficients, for every added line $(m,n)$, the coefficient $a_{mn}$ was set equal to the sum of $a_{\ell}^{\prime}s$ along the $m-n$ path, yielding $a_{2,5}=0.1936$, $a_{10,14}=0.0439$, $a_{7,12}=0.0419$. We kept the reference pressure at node $1$ and the compression ratios constant as in [12], and drew $1,500$ random gas injections $\mathbf{q}$. To construct these samples, we perturbed the benchmark injections $\mathbf{q}_{0}$ that lie in the range $[-15.61,22.01]$ by a standardized normal deviation. The injection at node $20$ was set to the negative sum of the remaining injections to get $\mathbf{1}^{\top}\mathbf{q}=0$ for all samples.

Using the modified meshed Belgian NGN of Fig. 1 and the random gas injections, we tested the exactness of (G2) and the performance of Algorithm 1. Not all of the random injections were feasible for the GF problem – some violated (4b) or (3b). Problem (G2) was infeasible for $876$ out of the $1,500$ random instances. Since (G2) is a relaxation of (G1), these instances are apparently infeasible for (G1) too. The performance of (G2) and Algorithm 1 was tested on the remaining $624$ gas injection instances. To evaluate the success of (G2) in solving (G1), we calculated the inexactness gap $G$ defined as

[TABLE]

for the pressures and flows obtained by (G2) and Algorithm 1.

The ranked inexactness gap for the feasible GF instances is shown by the first curve in Fig. 2. The gap was less than $10^{-3}$ for more than $97\%$ of the feasible instances, while the maximum gap over all instances was $0.009$. This corroborates that the proposed solver performs well even when Condition 2 is not met. Fig. 3 shows the running time for solving (G2) and Algorithm 1 over the $624$ feasible instances. The average (median) running time was 0.96 sec (0.89 sec).

Considering Condition 1, we used the fixed pressure at node $1$ and the pressures obtained at node $7$ for the feasible GF instances, we solved (G2) again. Although the hypothesis of Th. 3 does not hold anymore, the inexactness gap was found to be less than $10^{-3}$ for more than $94\%$ of the instances; see the second curve in Fig. 2. Thus, the tests reveal that the novel solver successfully finds the GF solution even when the sufficient Conditions 1–2 are violated. However, Condition 3 prohibiting gas circulations could not be violated for the Belgian NGN because the only active cycle in this NGN has parallel compressors, hence avoiding circulations from (4b). We next deal with GF instances on the GasLib-40 network, wherein a circulation could potentially occur.

GasLib-40 roughly represents a part of the German gas transport network [30]. The network exhibits 40 nodes, 39 pipes, and 6 compressors; see Fig. 1 (Bottom). The pipe dimensions, roughness coefficients, and a nominal demand vector $\mathbf{q}_{0}$ were derived from [30]. The goals for conducting additional tests on GasLib-40 include: i) Evaluating our solvers on a realistic setup; ii) Testing our MI-QCQP when Condition 3 is violated; and iii) Benchmarking the performance of our solvers against NR-based solver. We next briefly introduce the NR-based solver used for benchmarking. Given an injection $\mathbf{q}$, compressor ratios $\alpha_{\ell}$’s, and reference pressure $\psi_{1}$, stack the unknowns as ${\mathbf{y}=[\phi_{1},\dots,\phi_{L},\psi_{2},\dots,\psi_{N}]^{\top}}$. Define the equality constraints (2), (3a), and (4a) collectively as $\mathbf{g}(\mathbf{y})=0$. Given an initial estimate $\mathbf{y}_{0}$, the NR-based solver would iterate as

[TABLE]

where $t$ is the iteration count; matrix $\mathbf{J}(\mathbf{y}_{t})$ is the Jacobian of $\mathbf{g}(\mathbf{y})$ evaluated at $\mathbf{y}_{t}$; and $\mu$ a step size. A solution $\mathbf{y}^{\star}$ obtained on convergence of NR updates would be deemed feasible if the inequalities (3b) and (4b) are satisfied. Since, the NR updates target at attaining $\mathbf{g}(\mathbf{y})=\mathbf{0}$, the performance evaluation criteria for our results would be $\|\mathbf{g}(\mathbf{y})\|_{2}$ in lieu of the inexactness gap $G$.

In the first set of tests on GasLib-40, we generated $500$ gas injection instances $\mathbf{q}$ by scaling the entries of $\mathbf{q}_{0}$ independently, by random factors chosen uniformly on $[0.75,1.25]$. The pressure at node $1$ was set to $50$ bar and its injection was set to the negative sum of other nodes for all instances. Next, the compression ratios for the 6 compressors were drawn uniformly within $[1,2]$. All $500$ instances were solved using three approaches: a1) the MI-QCQP and Algorithm 1; a2) NR with flows initialized at $(\mathbf{A}^{\top})^{\dagger}\mathbf{q}$, and all pressures initialized at $\psi_{1}$; and a3) NR with flows and pressures initialized at the solution of MI-QCQP and Algorithm 1. The stopping criteria for NR was set to $\|\mathbf{g}(\mathbf{y})\|_{2}<10^{-3}$, subject to a maximum iteration count of $50$. The step size for both initialization scenarios was kept as $\mu=1$. The MI-QCQP deemed $5$ out of the $500$ instances as infeasible and the performance criteria $\|\mathbf{g}(\mathbf{y})\|_{2}$ was found to lie in $[0.005,0.183]$ with the median at $0.009$. To compare to the index of inexactness gap, the range for $G$ for the $495$ feasible cases was $[8\cdot 10^{-5},~{}6\cdot 10^{-2}]$. Thus, the MI-QCQP alongside Algorithm 1 was successful in finding the GF solution for all $495$ instances. Interestingly, $474$ of the $495$ feasible GF instances exhibit circulations, and hence violate Condition 3. Thus, the numerical results empirically demonstrate that the developed MI-QCQP alongside Algorithm 1 successfully solves the GF problem even when the conditions of Theorem 3 are violated. The NR solver, if initialized at the solution of MI-QCQP improves the solution accuracy, resulting in $\|\mathbf{g}(\mathbf{y}^{\star})\|_{2}$ within $1.2\cdot 10^{-4}-0.13$. For the $5$ instances deemed infeasible by MI-QCQP, the NR solver was initialized at all zero flows and pressures; all $5$ instances failed to converge. Surprisingly, when the NR solver was initialized with $(\mathbf{A}^{\top})^{\dagger}\mathbf{q}$ as flows and $\psi_{0}^{\prime}s$ as pressures, all $500$ instances failed to converge. The non-convergence of the NR solver is however alleviated when $\mu$ was reduced as discussed next.

A second set of tests were conducted on the GasLib-40 NGN with $500$ random injections and compressor ratios generated as described earlier. The MI-QCQP solver deemed $7$ of the $500$ instances as infeasible. All $500$ instances were then solved with the NR-based solver with flows initialized at $(\mathbf{A}^{\top})^{\dagger}\mathbf{q}$ and pressures at $\psi_{1}$, and $\mu$ was set to $\mu=0.9$. A steep decline in $\|\mathbf{g}(\mathbf{y}_{t})\|_{2}$ was observed in the first few (roughly 10) iterations, while the tolerance of $10^{-3}$ was not attained within the $50$ iterations limit. However, if the NR solver is initialized at the solution of MI-QCQP and Algorithm 1, the convergence criteria of $10^{-3}$ was attained at an average of $7.8$ iterations. The values of $\|\mathbf{g}(\mathbf{y})\|_{2}$ attained by three solution techniques a1)–a3) are shown in Fig. 4. The results suggest that the accuracy of the MI-QCQP solver is better than that of a2), which is a prudent initialization. However, if the NR-based solver is warm-started with the solution of MI-QCQP, an order of magnitude improvement in accuracy is observed. On the computational front, the MI-QCQP solver alongside Algorithm 1 is efficient with median solving time of $1.52$ sec. However, as anticipated, the NR solvers have superior performance with median solving time of $0.17$ sec. Finally, inspecting the $7$ instances deemed infeasible by the MI-QCQP solver, the solution obtained by a2) indicates violation of (4b); demonstrating the merit of the proposed MI-QCQP towards certifying infeasibility of GF instances.

VIII Conclusions

Exploiting recent results from graph theory and convex relaxations, this work provides a fresh perspective on the steady-state GF problem. The uniqueness of the GF solution has been established in a generalized setting for arbitrary NGN topologies, multiplicative compressors and multiple fixed-pressure nodes. Granted that the GF solution is unique, constrained and unconstrained versions of convex energy function minimization-based GF solvers have been proposed. These solvers can efficiently solve any GF task instance with a single fixed-pressure node and networks with compressors not on cycles. To expand the scope, an MI-QCQP GF solver had been also proposed relying on a convex relaxation of the Weymouth equation. The relaxation has been shown to be exact under specific network conditions. Numerical tests reveal that the developed MI-QCQP solver succeeds in finding the unique GF solution even when the needed conditions are violated. The success of the MI-QCQP relaxation is attributed to a judiciously designed objective. The developed approach sets forth an analytical platform for ensuring exact relaxation. Evaluating the performance of the developed approach for various optimal gas flow tasks constitutes an interesting research direction.

Proof of Lemma 2.

For an edge $\ell_{i}\in\mathcal{P}_{mn}$, let us name the incident node closer to $m$ as $m_{i}$, and the other node as $m_{i+1}$, as shown in Fig. 5.

Let $\boldsymbol{\psi}$ and $\boldsymbol{\psi}^{\prime}$ be the pressure vectors corresponding to $\boldsymbol{\phi}$ and $\boldsymbol{\phi}^{\prime}$. Since pressures $\psi_{m}$ and $\psi_{n}$ are fixed, it follows $\psi_{m}=\psi_{m}^{\prime}$ and $\psi_{n}=\psi_{n}^{\prime}$. Proving by contradiction, suppose (5a) holds. If that is the case, first it will be shown that $\psi_{m_{i}}^{\prime}-\psi_{m_{i+1}}^{\prime}>\psi_{m_{i}}-\psi_{m_{i+1}}$ for every lossy pipe $\ell_{i}\in\mathcal{P}_{mn}$.

Suppose that $\operatorname{sign}(\boldsymbol{\phi}^{\prime}-\boldsymbol{\phi})\cdot\boldsymbol{\pi}^{mn}>\mathbf{0}$. Let us denote the RHS of (3a) by $w(\phi_{\ell})$. It is evident that $w(\phi_{\ell})$ is monotonically increasing in $\phi_{\ell}$. Hence, for any lossy pipe $\ell_{i}\in\mathcal{P}_{mn}$, it holds

[TABLE]

where $(a)$ holds by hypothesis; $(b)$ stems from the monotonicity of $w(\phi_{\ell})$; $(c)$ holds because $\pi^{mn}_{\ell_{i}}\in\{0,1,-1\}$; $(d)$ holds from the property of $\operatorname{sign}$ by definition; and $(e)$ from the definition of $\boldsymbol{\pi}^{mn}$ and (3a). The inequality (VIII) implies

[TABLE]

Let us now apply (13) and (4a) for the edges $\ell_{1}$ to $\ell_{k}$ along $\mathcal{P}_{mn}$. For the fixed pressure node $m$, we have $\psi_{m}=\psi_{m}^{\prime}$. If $\ell_{1}$ is a lossy pipe, we get $\psi_{m_{2}}^{\prime}<\psi_{m_{2}}$ from (13); otherwise $\psi_{m_{2}}^{\prime}=\psi_{m_{2}}$ from (4a). Similarly, we can show that $\psi_{m_{3}}^{\prime}\leq\psi_{m_{3}}$, where the equality holds only if both $\ell_{1}$ and $\ell_{2}$ are compressors. However, this is practically impossible as every compressor is modeled as an ideal compressor followed by a lossy pipe, necessitating $\psi_{m_{3}}^{\prime}<\psi_{m_{3}}$. Continuing the process for all edges along $\mathcal{P}_{mn}$ yields $\psi_{n}^{\prime}<\psi_{n}$, which contradicts with node $n$ being a fixed-pressure node. Similarly, the assumption $\operatorname{sign}(\boldsymbol{\phi}^{\prime}-\boldsymbol{\phi})\cdot\boldsymbol{\pi}^{mn}<0$ leads to a contradiction by yielding $\psi_{n}^{\prime}>\psi_{n}$. ∎

Proof of Lemma 3.

Given the two pairs $(\mathbf{q},\boldsymbol{\phi})$ and $(\mathbf{q}^{\prime},\boldsymbol{\phi}^{\prime})$ satisfying (2) and $\mathbf{q}\neq\mathbf{q}^{\prime}$, let us define $\tilde{\boldsymbol{\phi}}:=\boldsymbol{\phi}^{\prime}-\boldsymbol{\phi}$ and $\tilde{\mathbf{q}}:=\mathbf{q}^{\prime}-\mathbf{q}$. By applying (2) on $(\mathbf{q},\boldsymbol{\phi})$ and $(\mathbf{q}^{\prime},\boldsymbol{\phi}^{\prime})$, and taking the difference, we get

[TABLE]

Since $\mathbf{1}\in\operatorname{null}(\mathbf{A})$, premultiplying (14) by $\mathbf{1}^{\top}$ provides

[TABLE]

From (14)–(15), the pair $(\tilde{\mathbf{q}},\tilde{\boldsymbol{\phi}})$ qualifies as a set of balanced gas injections. By definition of $(\tilde{\mathbf{q}},\tilde{\boldsymbol{\phi}})$, proving (6) is equivalent to showing there exists a path $\mathcal{P}_{mn}$ for which

[TABLE]

To prove the existence of such a path, we use the ensuing result based on [31, Th. 8.8].

Lemma 5** ([31]).**

Given a graph with injection $q$ at node $s$, demand $q$ at node $t$, and zero injections at all other nodes, there exists an $s$-$t$ path with flow directions along the path from $s$ to $t$.

Lemma 5 considers a single-source single-destination network flow setup. We transform our problem to this setup through the next steps; see also Fig. 6:

1. 1)

The nodes of graph $\mathcal{G}$ are partitioned into the subset with positive $\mathcal{N}_{+}:\{n\in\mathcal{N}:\tilde{q}_{n}>0\}$; negative $\mathcal{N}_{-}:\{n\in\mathcal{N}:\tilde{q}_{n}<0\}$; and zero injections $\mathcal{N}_{0}:\{n\in\mathcal{N}:\tilde{q}_{n}=0\}$. Because $\tilde{\mathbf{q}}\neq\mathbf{0}$, the sets $\mathcal{N}_{+}$ and $\mathcal{N}_{-}$ are non-empty. 2. 2)

Augment $\mathcal{G}$ by adding nodes $s$ and $t$. 3. 3)

All nodes in $\mathcal{N}_{+}$ are connected to node $s$, and all nodes in $\mathcal{N}_{-}$ are connected to node $t$. 4. 4)

The injections in $\mathcal{N}_{+}$ are lumped in node $s$ by setting the flows $\tilde{\phi}_{sn}=\tilde{q}_{n}$ for all $n\in\mathcal{N}_{+}$. Similarly, the demands in $\mathcal{N}_{-}$ are lumped in node $t$ by setting the flows $\tilde{\phi}_{nt}=-\tilde{q}_{n}$ for all $n\in\mathcal{N}_{-}$.

Applying Lemma 5 on this augmented graph, there exists a path $\mathcal{P}_{st}$ with flow directions from $s$ to $t$. For any such path $\mathcal{P}_{st}$, eliminate the first and last edges to get a path $\mathcal{P}_{mn}$ with $m\in\mathcal{N}_{+}$ and $n\in\mathcal{N}_{-}$. Claim (16b) follows by construction. We next show (16a): For each edge $\ell\in\mathcal{P}_{mn}$, it was shown that the direction of $\tilde{\phi}_{\ell}$ is along the path $\mathcal{P}_{mn}$. If $\pi_{\ell}^{mn}=+1$, the direction of edge $\ell$ agrees with the direction of $\mathcal{P}_{mn}$. Since $\tilde{\phi}_{\ell}$ is along $\mathcal{P}_{mn}$, then $\tilde{\phi}_{\ell}>0$. If $\pi_{\ell}^{mn}=-1$, the direction of edge $\ell$ is opposite to the direction of $\mathcal{P}_{mn}$. Since $\tilde{\phi}_{\ell}$ is along $\mathcal{P}_{mn}$, then $\tilde{\phi}_{\ell}<0$. Either way, it holds that $\tilde{\phi}_{\ell}\pi_{\ell}^{mn}>0$ for all $\ell\in\mathcal{P}_{mn}$, which proves (6a). ∎

Proof of Theorem 3.

Before proving the main result, we will need two preliminary results.

Lemma 6**.**

For a lossy pipe $\ell=(m,n)$ not on an active cycle, if the triplet $(\psi_{m},\psi_{n},\phi_{\ell})$ satisfies (11), then the triplet $(\psi_{m}+\delta,\psi_{n}+\delta,\phi_{\ell})$ also satisfies (11) for any finite $\delta$.

Lemma 6 follows directly from the fact that (11) involves pressure differences rather than pressures.

Lemma 7**.**

Consider an active cycle $\mathcal{C}_{0}$ and index its nodes as $\{0,\dots,k\}$. Given a fixed pressure $\psi_{0}$ and flows $\{\phi_{\ell}\}_{\ell\in\mathcal{C}_{0}}$ satisfying Condition 3 and (4b), there exists a set of pressures $\{\psi_{i}\}_{i=1}^{k}$ satisfying (11) and (4a).

Proof.

From Condition 3 and the fact that a compressor is modeled as an ideal compressor followed by a lossy pipe, it is not hard to see that there must exist a node $k\in\mathcal{C}_{0}$ that leads to one of the four flow scenarios shown in Fig. 7.

Proving by construction, we will next define pressures $\{\psi_{i}\}_{i=1}^{k}$ such that (11) and (4a) are satisfied for all edges in $\mathcal{C}_{0}$. Traversing the paths $0\rightarrow k_{1}$ and $0\rightarrow k_{2}$, one can recursively define pressures for all nodes using $\psi_{0}$ and flows $\{\phi_{\ell}\}_{\ell\in\mathcal{C}_{0}}$ based on the exact Weymouth equation (3) and (4a). The pressures on the remaining nodes of $\mathcal{C}_{0}$ can be defined for the four scenarios of Fig. 7 as follows:

[TABLE]

To see that the constructed pressures satisfy (11), take for example scenario $(a)$. Applying (11) along the edges $(k_{1},k)$ and $(k_{2},k)$ yield that $\psi_{k}$ should satisfy $\psi_{k}\leq\psi_{k_{1}}-a_{k_{1}k}\phi_{k_{1}k}^{2}$ and $\psi_{k}\leq\psi_{k_{2}}-a_{k_{2}k}\phi_{k_{2}k}^{2}$. This is indeed the case by selecting $\psi_{k}$ as the minimum of the two RHS. Similar reasoning applies to the other scenarios. ∎

Proceeding with the proof of Theorem 3, let $(\boldsymbol{\phi},\boldsymbol{\psi})$ be the unique solution to (G1), and $(\boldsymbol{\phi}^{\prime},\boldsymbol{\psi}^{\prime})$ a minimizer of (G2). Proving by contradiction, assume that there exists an edge $\ell$ not belonging to an active cycle, such that $\phi_{\ell}^{\prime}\neq\phi_{\ell}$. Recall that the set of all active cycles is denoted by $\mathcal{S}_{\mathcal{C}}^{a}$. Since both flow vectors satisfy (2), their difference $\mathbf{n}:=\boldsymbol{\phi}-\boldsymbol{\phi}^{\prime}$ must lie in the nullspace of $\mathbf{A}^{\top}$. The nullspace of $\mathbf{A}^{\top}$ is spanned by the indicator vectors for all fundamental cycles in the gas network graph [32, Corollary 14.2.3]. Therefore, the entries of $\mathbf{n}$ related to edges not on a cycle must be zero. Since by hypothesis $\ell\notin\mathcal{S}_{\mathcal{C}}^{a}$, edge $\ell$ should belong to one of the cycles in $\mathcal{S}_{\mathcal{C}}\setminus\mathcal{S}_{\mathcal{C}}^{a}$. This non-active cycle will be henceforth termed $\mathcal{C}$.

The rest of the proof is organized in three parts: Part I constructs a flow vector $\hat{\boldsymbol{\phi}}$ that satisfies (2) and (4b). Part II shows there exists a $\hat{\boldsymbol{\psi}}$ so that the pair $(\hat{\boldsymbol{\phi}},\hat{\boldsymbol{\psi}})$ is feasible for (G2). Part III shows that $(\hat{\boldsymbol{\phi}},\hat{\boldsymbol{\psi}})$ attains a smaller objective for (G2), thus contradicting the optimality of $(\boldsymbol{\phi}^{\prime},\boldsymbol{\psi}^{\prime})$.

Part I: Define the flow vector $\hat{\boldsymbol{\phi}}$ as

[TABLE]

By construction, vector $\hat{\boldsymbol{\phi}}$ satisfies

[TABLE]

where $\mathbf{n}^{\mathcal{C}}$ is the indicator vector for cycle $\mathcal{C}$; the constant $\lambda$ is nonzero; and vector $\mathbf{n}^{a}\in\operatorname{null}(\mathbf{A}^{\top})$ can have nonzero entries only for edges in active cycles. Since $\boldsymbol{\phi}^{\prime}$ satisfies constraint (2) and $\mathbf{A}^{\top}\mathbf{n}^{\mathcal{C}}=\mathbf{A}^{\top}\mathbf{n}^{a}=\mathbf{0}$, then $\mathbf{A}^{\top}\hat{\boldsymbol{\phi}}=\mathbf{A}^{\top}\boldsymbol{\phi}^{\prime}=\mathbf{q}$. This proves that $\hat{\boldsymbol{\phi}}$ satisfies (2). Note that $\hat{\boldsymbol{\phi}}$ is constructed by selecting entries from $\boldsymbol{\phi}$ and $\boldsymbol{\phi}^{\prime}$. Granted both $\boldsymbol{\phi}$ and $\boldsymbol{\phi}^{\prime}$ satisfy (4b), vector $\hat{\boldsymbol{\phi}}$ trivially satisfies (4b) too.

Part II: We will delineate the steps for constructing a vector of pressures $\hat{\boldsymbol{\psi}}$ such that $(\hat{\boldsymbol{\phi}},\hat{\boldsymbol{\psi}})$ is feasible for (G2). Let us select a spanning tree $\mathcal{T}$ of the NGN graph $\mathcal{G}$ rooted at the reference $r$. We shall define the pressures $\hat{\psi}_{n}$’s while traversing $\mathcal{T}$ via depth-first search. In such a traversal, the following three cases may be identified on arriving at any node $n$:

Case 1: Node $n$ is neither in $\mathcal{C}$ nor on an active cycle. Let $n-1$ be the parent node of $n$ in $\mathcal{T}$ and define

[TABLE]

Since the edge $(n-1,n)$ is not in $\mathcal{C}\cup\mathcal{S}_{\mathcal{C}}^{a}$, we have $\hat{\phi}_{n-1,n}={\phi}_{n-1,n}^{\prime}$ from (17). Therefore, if $(n-1,n)$ is a lossy pipe, Lemma 6 ensures that the defined pressure $\hat{\psi}_{n}$ satisfies (11). Moreover, if $(n-1,n)$ is a compressor, constraint (4a) is satisfied trivially by definition.

Case 2: Node $n$ is in $\mathcal{C}$. If $n$ is the first node in $\mathcal{C}$ to be visited, define $\hat{\psi}_{n}$ as in Case 1. Then, define the pressures for the remaining nodes $i\in\mathcal{C}$ as $\hat{\psi}_{i}:=\psi_{i}+(\hat{\psi}_{n}-\psi_{n})$. Note from (17) that the flows along $\mathcal{C}$ are assigned from $\boldsymbol{\phi}$, the pair $(\boldsymbol{\phi},\boldsymbol{\psi})$ satisfies (3) and hence the relaxed Weymouth (11) as well. The constructed pressures $\hat{\psi}_{i}$’s for $i\in\mathcal{C}$ are simply a shifted version of the pressures $\psi_{i}$’s. Therefore, the pressures $\hat{\psi}_{i}$’s satisfy (11) from Lemma 6. Mark all nodes in $\mathcal{C}$ as traversed and continue.

Case 3: Node $n$ is in an active cycle $\mathcal{C}_{a}$. If $n$ is the first node in $\mathcal{C}_{a}$ to be traversed, define the $\hat{\psi}_{n}$ as in Case 1. Then, define the pressure for the remaining nodes $i\in\mathcal{C}_{a}$ using Lemma 7. Mark all nodes in $\mathcal{C}_{a}$ as traversed and continue.

Since the constructed pressures satisfy (11) and (4a), the pair $(\hat{\boldsymbol{\phi}},\hat{\boldsymbol{\psi}})$ is feasible for (G2). Observe that the pressure drop across lossy pipes not in $\mathcal{C}$ is ${\hat{\psi}_{m}-\hat{\psi}_{n}=\psi_{m}^{\prime}-\psi_{n}^{\prime}}$ for Case 1; and ${\hat{\psi}_{m}-\hat{\psi}_{n}=\psi_{m}-\psi_{n}}$ for lossy pipes in $\mathcal{C}$ under Case 2. This fact is imperative for the ensuing Part III.

Part III: We will next show that $r(\boldsymbol{\psi}^{\prime})>r(\hat{\boldsymbol{\psi}})$ to contradict the optimality of $\boldsymbol{\psi}^{\prime}$. Note that the objective $r(\boldsymbol{\psi})$ in (G2) sums up the absolute pressure differences along lossy pipes, but not on active cycles. Since by construction these differences have changed only along $\mathcal{C}$, we get

[TABLE]

As the pressure differences depend on flows, we next compare the entries of $\hat{\boldsymbol{\phi}}$ and $\boldsymbol{\phi}^{\prime}$ along $\mathcal{C}$ using (18). Since the edge directions are assigned arbitrarily, assume wlog that $\hat{\phi}_{mn}\geq 0$ for all $(m,n)\in\mathcal{C}$. Given $\mathbf{n}_{\mathcal{C}}$ and (18), one can find the value of $\lambda$. If $\lambda<0$, reverse the reference direction for cycle $\mathcal{C}$ to get a positive $\lambda$. Because of this, we can assume $\lambda>0$.

Recall that $\mathbf{n}_{\mathcal{C}}\in\{0,\pm 1\}^{P}$. Partition the set of edges in $\mathcal{C}$ into mutually exclusive sets $\hat{\mathcal{P}}_{+}$ and $\hat{\mathcal{P}}_{-}$ based on positive and negative entries of $\mathbf{n}_{\mathcal{C}}$, respectively. From (18), it follows

[TABLE]

Summing up the pressure drops along $\mathcal{C}$ for $\hat{\boldsymbol{\psi}}$ should be zero. Since the pressure drops along $\mathcal{C}$ are positive for the edges in $\hat{\mathcal{P}}_{+}$, and negative along the edges in $\hat{\mathcal{P}}_{-}$, it holds that

[TABLE]

where the absolute value is trivial since $\hat{\phi}_{mn}\geq 0$ for all $(m,n)\in\mathcal{C}$.

Drawing similar relations on $\boldsymbol{\psi}^{\prime}$, define the set $\mathcal{P}_{+}^{\prime}\subset\mathcal{C}$ containing any edge $(m,n)\in\mathcal{C}$ such that the flow $\phi_{mn}^{\prime}$ is along the direction of $\mathbf{n}_{\mathcal{C}}$. Using the same argument as in (VIII) for $\boldsymbol{\psi}^{\prime}$, we obtain

[TABLE]

Because the flows in $\hat{\boldsymbol{\phi}}$ for the edges in $\hat{\mathcal{P}}_{+}$ are aligned with $\mathbf{n}_{\mathcal{C}}$ and $\phi_{\ell}^{\prime}>\hat{\phi}_{\ell}$ for these edges from (20), it follows that $\hat{\mathcal{P}}_{+}\subseteq\mathcal{P}_{+}^{\prime}$. Using the latter in (22), we get

[TABLE]

For every edge $\ell=(m,n)\in\hat{\mathcal{P}}_{+}$, it holds that

[TABLE]

where $(a)$ comes from the definition of pressures in Case 2 of Part II; $(b)$ descends from $\phi_{\ell}^{\prime}>\hat{\phi}_{\ell}>0$; and $(c)$ from (11). Summing (24) over all $\ell\in\hat{\mathcal{P}}_{+}$ and multiplying by 2 gives

[TABLE]

where the inequality stems from (22) and (VIII). From (19), the latter implies that $r(\boldsymbol{\psi}^{\prime})>r(\hat{\boldsymbol{\psi}})$, hence contradicting the optimality of $\boldsymbol{\psi}^{\prime}$. ∎