Combinatorial Acyclicity Models for Potential-based Flows

Potential-based ﬂows constitute a basic model to represent physical behavior in networks. Under natural assumptions, the ﬂow in such networks must be acyclic. The goal of this paper is to exploit this property for the solution of corresponding optimization problems. To this end, we introduce several combinatorial models for acyclic ﬂows, based on binary variables for ﬂow directions. We compare these models and introduce a particular model that tries to capture acyclicity together with the supply/demand behavior. We analyze properties of this model, including variable ﬁxing rules. Our computational results show that the usage of the corresponding constraints speeds up solution times by about a factor of 3 on average and a speed-up of a factor of almost 5 for the time to prove optimality.


Introduction
Potential-based flows form a basic model for physical networks.Such flows are necessarily acyclic, which we will exploit in this article.To introduce the basic idea, let D = (V, A) be a simple directed graph without anti-parallel arcs.For all arcs a ∈ A, there is a continuous, strictly increasing potential function ψ a : R R with ψ a (0) = 0 as well as a resistance β a > 0. Each node v ∈ V has an associated potential π v .Under mild assumptions on the potential functions and given demand on the nodes, the defining equations induce a unique flow x ∈ R A and unique potential differences (see Section 2).For physical networks, the potentials π u correspond to quantities like squared pressure or voltage.Note that a directed instead of undirected graph D is used to define a direction of the flow and (1).Thus, flow values can also be negative, indicating flow in the opposite direction of the arc.
A flow x ∈ R A defines a directed graph D(x) = (V, A(x)) with If x satisfies (1), this graph is always acyclic.To see this, assume that there would exist a directed cycle C ⊆ A(x).Then splitting the arcs in C into forward arcs, i.e., those contained in the original graph, and backward arcs, i.e., those which have the opposite direction to the original graph, we obtain where we use that C is a cycle and thus the alternating sum of the potentials vanishes.However, since the resistances β a are positive and each ψ a is strictly increasing with ψ a (0) = 0, we have β a ψ a (x a ) > 0 if x a > 0 and β a ψ a (x a ) < 0 if x a < 0. Thus, the value of the first line is positive, leading to a contradiction.This shows that D(x) is acyclic, corresponding to the physical property of having a conservative potential.Uniqueness and acyclicity are two important physical properties that are captured by potential-based flows.Moreover, this model class is important for handling energy networks; see Section 2 for examples.Such networks often contain active network elements like switches/valves, allowing to close a connection, or generators/pumps/compressors, which can increase the potential on certain arcs.These elements allow to control flow and potentials.Their presence may violate acyclicity, which provides more freedom to control the network.However, the passive components remaining after removal of active elements satisfy acyclicity.
Using the degrees of freedom of active elements, several different optimization problems over such networks are interesting, e.g., energy minimal operation under the assumption that a flow demand is satisfied.When solving such optimization problems to global optimality, one can exploit the fact that the passive components of the network still have an acyclic flow.This is the main idea of this article.We will demonstrate how to enhance existing mixed-integer nonlinear programing (MINLP) formulations using binary variables for the flow directions and constraints that enforce acyclicity.To this end, we introduce a nested sequence of polytopes that encode the directions of the flows, relaxing more and more constraints of the nonlinear model of potential-based flows along the way.Each of these polytopes provides a combinatorial model of acyclic flows.We investigate a particular model that provides a good compromise between the nonlinear model and the so-called acyclic subgraph polytope in more detail.We then add the corresponding inequalities to a potential-based network optimization problem for gas network examples.We demonstrate that this approach leads to an improvement of the solving time by about a factor of 3 on average and a significant speed-up for the time to prove optimality by almost a factor of 5.
This article is structured as follows.In Section 2 we discuss potential-based flows in more detail and provide examples.Section 3 first introduces the sequence of combinatorial models.Their relation is studied in Section 3.1.In Section 3.2 we investigate a model solely based on acyclic directions and the corresponding computational complexity.Section 3.3 then introduces our main model, exploiting both acyclicity and the fact that one needs to connect sources and sinks.This model is investigated in Sections 3.4 and 3.5 in more detail.Then Section 4 presents the computational results.We close with a conclusion and some open questions in Section 5.
Literature Review Potential-based flows have been used in many different contexts.Hendrickson and Janson [18] provide an overview.We will refer to more literature in Section 2 and also discuss examples.The special case of gas transport will be used in our computations in Section 4.
The topic of acyclic flows for potential-based flow has been investigated by Becker and Hiller in four articles [19,2,3,4].Their motivation is similar to ours and they also test their methods on gas networks.The main combinatorial model of these articles is based on socalled acyclic source transshipment sink (ASTS) orientations.We will arrive at an equivalent definition through a polyhedral approach in Section 3.3.Their contributions can be briefly summarized as follows.A characterization of the cases in which an ASTS orientation exists is given in [3].Moreover, various decomposition results based on 2-connected components are presented in [2,4].This allows to preprocess the networks [2, Section 3].The enumeration of ASTS orientations is discussed in [2,4].This is used in a Dantzig-Wolfe type approach to strengthen potential-based flow formulations for gas networks [2,4].
Our results differ from the ones by Becker and Hiller in the following way.We embed the common combinatorial model in a sequence of polytopes, which we investigate polyhedrally.We use a different binary encoding, provide a different setup, investigate complexity results, and add the inequalities to the model instead of using enumeration of the configurations.The ideas of this article were used for the computations of [17] and we provide the background in this article.Apart from the mentioned literature, we are not aware of any other works concerning combinatorial models for acyclic flows.

Potential-based Flows
Potential-based flows have been studied repeatedly in the literature.It seems that the first appearance is in Birkhoff and Diaz [6].A general treatment appears in Rockafellar [33].Before referring to more results in the literature, we first complete the setting of potentialbased flows.
For every node v ∈ V, there are lower and upper pressure bounds π v and π v ∈ R, respectively, with π v ≤ π v .Additionally, for every arc a ∈ A, there are lower and upper flow bounds x a and x a ∈ R, respectively, with x a ≤ x a .Let b ∈ R V be a supply and demand vector that is balanced, i.e., v∈V b satisfies the the following constraints: We call x ∈ R A a b-flow if it satisfies (2a).Throughout this paper, we will assume that each potential function ψ a : R R, a ∈ A, is continuous and strictly increasing with ψ a (0) = 0.For some results in the literature, additional requirements on the potential functions are needed, e.g., that they are odd (i.e., ψ a (x) = −ψ a (−x)), positively homogeneous (i.e., ψ a (λx) = λ r ψ a (x) for all x ∈ R and λ > 0 with some constant r > 0) or that they are the same for every arc, see, e.g., [15].
One important result, see Maugis [27], Collins et al. [8], and Ríos-Mercado et al. [32] is the following: Assume that D is weakly connected, there are no bounds on the potentials and flows, for a given node s ∈ V the potential π s is fixed, and the potential functions are continuous and strictly increasing.Then there exists a unique feasible potential-based flow (x, π).Consequently, System (2) with a fixed potential π s is either infeasible or has a unique solution.One tool for proving uniqueness is a cost minimal flow problem with a strongly convex objective, whose dual multipliers provide the potentials, see [27], [33], and [15] for more information and a discussion of the corresponding Lagrange dual.
Example 1. Several interesting applications can be modeled as potential-based flows, see, e.g., Hendrickson and Janson [18].We present three energy network examples: 1. Stationary Gas Transport Networks: Here arcs correspond to pipe(lines), the potentials are the squares of pressures, and flows are gas mass flows.One common approximation of gas flow is (1) with ψ a (x a ) = |x a | x a .The positive arc-specific constant β a depends on the pipe diameter, length, and roughness of its inner wall.More details on stationary gas flow in pipeline networks are given in the book [22] and the references therein.
The above model assumes constant heights of the network, but one can use scaling to incorporate different heights, see [15].
The computational results in Section 4 are based on gas networks, extended by active elements like valves and compressors.
2. Water Networks: Here potentials correspond to hydraulic heads.3. Lossless DC Power Flow Networks: In this case, the potentials are voltages and the potential function is linear ψ a (x a ) = x a .For more information about power flow network planning, we refer to [5].

Combinatorial Models for Acyclic Flows
To study acyclic flows, we first introduce some basic notation.Consider the simple and weakly connected directed graph D without anti-parallel arcs.Note that these assumptions are without loss of generality: Loops always have a zero flow and can be removed, anti-parallel arcs can be reoriented to be parallel and parallel arcs can be merged into one arc with adapted β-value.Moreover, each weakly connected component can be treated separately.We will use a to denote the reverse arc of some arc a ∈ A, i.e., if a = (u, v), then a := (v, u).Furthermore, let A := { a : a ∈ A} be the set of the reversed arcs.Note that for a given flow x, we can also write and the digraph D(x) = (V, A(x)) is a reorientation of the subgraph consisting of arcs with nonzero flow such that all arcs point in the direction of the flow.We say that x is an acyclic flow if D(x) is acyclic.An alternative definition for acyclic flows is given by Hiller and Becker [19,Definition 1].
Remark 2. Note that being an acyclic s-t-flow is not related to having a flow decomposition which only consists of paths and no cycles.On the one hand, the sum of flow along paths can contain a cycle.For example, consider the flow network given in Figure 1.The s-t-flow depicted in this figure can be decomposed into the flows x 1 (dashed) and x 2 (dotted).Both x 1 and x 2 are flows along paths and therefore acyclic, nevertheless their sum is not, since D(x 1 + x 2 ) contains the directed cycle C = (u, v, w, u).On the other hand, the sum of flows along paths and cycles can be acyclic.For example consider flow x 3 with value −1 along C.Then, D(x 1 + x 2 + x 3 ) is acyclic.
The following example shows how flow directions can depend on the resistances β a .Example 3. Consider the potential network given in Figure 2 with source s and sink t and By constraint (2a), it is clear that x 1 + x 2 = x 4 + x 5 = b s > 0 has to hold for every potential flow (x, π).Furthermore, x 1 , x 2 , x 4 and x 5 have to be nonnegative: Assume that one of them is negative.We can assume w.l.o.g.x 2 < 0 by symmetry.Then x 1 > 0 by flow conservation.Having x 3 > 0 would close the cycle (s, u, v, s).Thus, x 3 ≤ 0 and by flow conservation x 5 < 0. Since also x 4 > 0, this closes a cycle (s, u, t, v, s).
Thus, except for arc (u, v), all flow directions for this potential network are fixed, independent of the ψ a and the β-values.However, the flow direction of (u, v) depends on the ψ a and the β-values.Indeed, if all β-values and all ψ = ψ a are equal, then the corresponding flow is x 3 = 0.Moreover, consider the case which is a contradiction.Hence, x 3 < 0. Furthermore, by symmetry To express acyclicity, for each arc a ∈ A, we introduce binary variables z + a and z − a that model the flow direction as follows: where sgn(x a ) = −1 if x x < 0, sgn(x a ) = 0 if x a = 0, and sgn(x a ) = 1 if x a > 0. Thus, these constraints imply that z + a = 1 if x a > 0, z − a = 1 if x a < 0, and z + a = z − a = 0 if x a = 0.The total model for potential-based flows is the following: and the corresponding subgraph In order to exploit acyclicity of potential-based flows, we will investigate purely combinatorial models of acyclicity, i.e., polytopes that are only based on the variables z + a and z − a .There are several possibilities for such models, depending on how many properties of the potential-based flow are used.We present four polytopes in the following and one model in Section 3.3.The main goal is to derive inequalities that can be added to (PBF) in order to improve the computational solving performance.
Projected Potential-Based Flows The most specific model considers the projection of feasible points of (PBF) for a given network with given balanced b ∈ R V and yields the polytope of potential-based flow directions Note that since the flows are unique, (z + , z − ) is also unique.The only possible variation is whether P PF is empty or not.Since it is an NP-hard problem to decide whether there exists a potential-based flow for the case of DC-flows, see Lehmann et al. [26], and for the case of gas networks, see Szabó [38], it is an NP-hard problem to decide whether P PF is empty.
Projected Universal Potential-Based Flows Example 3 shows that the flow directions can depend on the values of β.We therefore investigate a model in which the resistances β are allowed to vary.Consider the asymptotic behavior of For fixed potentials and β a ∞ we get x a 0. Thus, we identify (2b) for β a = ∞ with the constraint x a = 0 and decoupled potentials π u , π v .This has the same effect as if arc a = (u, v) would not exist.In the following we denote the extended real line with Again given a digraph D with balanced b ∈ R V , we define the polytope of universal potential-based flow directions for (PBF) as By allowing the resistances β to vary, P UPF abstracts from the particular network to some extent.Note that the polytope is universal in the sense that changes in β allow a corresponding change of direction of some arcs as in Example 3.

Acyclic Flows
The polytope of acyclic flow directions is Note that the potential equation (2b) of (PBF) is relaxed.Because of (3), we could replace the requirement that D(x) is acyclic by acyclicity of D(z + , z − ).

Acyclic Subgraphs
The polytope that abstracts the most from potential-based flows is the polytope of acyclic subgraphs Note that acyclicity of D(z + , z − ) implies z + a + z − a ≤ 1 for each arc a ∈ A. In Section 3.3, we will refine this model by using knowledge of sources, sinks, and inner nodes.
Remark 4.There are two alternatives to using (3).The first one is These constraints form a relaxation of (3), since the direction variables z + a and z − a can be chosen freely if x a = 0 (as long as z + a + z − a ≤ 1).For (x, π, z + , z − ) ∈ X this would allow D(z + , z − ) to have cycles although D(x) is acyclic.Thus, P PF and P UPF would include (z + , z − ) that do not correspond to potential-based flows.Note that this model is only valid if x ≤ 0 ≤ x, because positive bounds x or negative bounds x would be overruled when setting z + = 0 or z − = 0.
The second alternative for (3) is This would require to assign a direction to a 0-flow arc, which and would make the following analysis more difficult.

Relations among Combinatorial Models
We begin with the obvious observation that In fact, Example 3 shows that in general the first inclusion is strict.The last inclusion is always strict: on the one hand, if b = 0, then 0 ∈ P AF , but we always have 0 ∈ P AS , on the other hand, if b = 0 then P AF = {0}, but we can always set a single direction to 1 in P AS .Moreover, we have: If there are no potential bounds then P UPF = P AF holds.
Proof.If suffices to show P AF ⊆ P UPF .To this end, suppose that b = 0. Then the only acyclic flow satisfying (2a) is x = 0. Thus, if x = 0 satisfies the flow bounds (2d), z + = z − = 0 is the only possible solution and P UPF = P AF = {0} holds.Otherwise, let b = 0 be balanced and consider any acyclic x ∈ R A that satisfies flow conservation (2a) and the bounds (2d).We first possibly reorient the graph and invert x such that x ≥ 0. Define z + a = 1 iff x a > 0 and z − a = 0 for all arcs a.Thus, (x, z + , z − ) satisfies (3) and we want to show that (z + , z − ) ∈ P UPF , i.e., there exists potentials π v that satisfy (2b) for appropriately chosen β a .
We first compute potentials π v such that π u − π v ≥ x a for all arcs a = (u, v) ∈ A and π u − π v = 0 if x a = 0.This can be done by adding an artificial node r and arcs (r, v) with weight 0 for all v ∈ V.The weights on the other arcs are given by −x.The graph does not contain negative cycles, because x is acyclic.Thus, the Moore-Bellman-Ford algorithm computes the shortest distances π v from r to v ∈ V with π v ≤ π u − x a for all a = (u, v) ∈ A, see, e.g., Korte and Vygen [23].Note that one can shift the potentials such that π ≥ 0 holds.
We now choose β-values such that This yields the following result.
Corollary 6.If no potential bounds are present, P PF ⊆ P UPF = P AF ⊂ P AS and the two inclusions are strict in general.
The following results justify the choice of β a = ∞ in P UPF .Therefore, given a weakly connected digraph D = (V, A) and a balanced supply and demand vector b ∈ R V , we define the set of all potential-based flows We first show that in the absence of flow and potential bounds the closure of X x is given by permitting β a = ∞.Note that X x is never empty.
Proposition 7. Let the potential function ψ be continuous and strictly increasing.Consider a weakly connected digraph D = (V, A) and a balanced supply and demand vector b ∈ R V .The closure of potential-based flows is given by Proof.In the case b = 0, both sets only contain x = 0. Thus, it suffices to consider the case b = 0.In the following, we denote the set on the right of ( 6) by X ∞ .We first show cl(X x ) ⊆ X ∞ .We assume there exists x * ∈ cl(X x ) \ X x as otherwise there is nothing to show.To prove x * ∈ X ∞ , we have to construct resistances β ∈ RA >0 and potentials π * such that (x * , π * ) satisfy (2a) and (2b).If x * is acyclic, we can use the procedure in the proof of Lemma 5. Thus, we have to show that A(x * ) is acyclic.
Suppose that A(x * ) contains a directed cycle C and w.l.o.g.assume that x * a > 0 for all a ∈ C. Let (x k ) k∈N ⊂ X x converge to x * .Then for some k 0 , x k 0 a > 0 holds for all a ∈ C, which contradicts x k 0 being a potential-based flow.Hence, x * is acyclic and we conclude that the inclusion cl(X x ) ⊆ X ∞ holds.
We now show the reverse inclusion , then we are done due to definition.Therefore, assume that β a = ∞ for at least one arc a ∈ A. We have to construct a sequence (x k ) k∈N ⊂ X x converging to x * to finish the proof.We construct this sequence as follows: 1. We choose any sequence (α k ) k∈N ⊂ R >0 with α k ∞.
2. We define a sequence of resistances (β k ) k∈N ⊂ R A >0 , by using β k a = β * a for all arcs a with β * a < ∞ and β k a = α k otherwise.3. Choose a source s ∈ V and fix π k s = π * s for all k ∈ N. 4. Since there are no flow or potential bounds, there exists a unique solution (x k , π k ) for every k ∈ N.
We claim that the sequence x k constructed this way converges to x * .To see this, consider the subgraph which results from removing all arcs with β * a = ∞ and all resulting isolated nodes.
Due to flow conservation and the acyclicity of potential-based flows, we can derive lower and upper bounds on the flow of each arc, e.g., the flows are bounded by the sum of inflows.These flow bounds together with (2b) then define a lower and an upper bound on the potential difference between the nodes of each weakly connected component in D ∞ .
Let a = (u, v) ∈ A\A ∞ and suppose that x k a does not converge to 0. Then (a subsequence of) the sequence β k a ψ a (x k a ) converges to ±∞.If u, v are in the same connected component of D ∞ , this contradicts the potential differences of each component being bounded.Otherwise, note that due to flow conservation there can only be flow from one connected component to another, if there is also flow to the first component from another component, and vice versa.Thus, there exists a "cycle" of arcs in A \ A ∞ connecting different components of D ∞ .Hence, if the flow on these arcs does not converge to 0, the potential differences of the nodes where the flow enters and leaves the different components converges to ±∞, which is a contradiction as before.Therefore, x k x * holds, which concludes the proof.
The previous result can also be extended to the case with flow bounds.
Corollary 8. Using the assumptions of Proposition 7, the following holds.Given flow bounds x ≤ x, we define X Proof.The inclusion "⊆" holds by the same arguments as before.To see "⊇", we use the same construction to define the sequence (x k , π k ) k∈N .After possibly choosing a subsequence, all elements of the sequence either satisfy the flow bounds, or all violate the flow bounds.That is, the sequence is either contained in cl(X ).In the first case, we are done.Otherwise, note that the flow bounds are satisfied in the limit and thus the limit is not contained in the complement of cl(X [x] ) but in the intersection of the closures.
• Combining Lemma 5 and Corollary 8 yields that in the absence of potential bounds, acyclic flows coincide with the closure of potential-based flows.
• Instead of taking the closure of the flows only, we could also consider Then the closure satisfies where additionally β a = 0 is permitted.Here, the reverse inclusion is in general not true, because when using β a = 0, flow in a cycle is possible, while potential-based flows are always acyclic.
• Note that taking the closure of potential-based flows together with the corresponding directions defined by (3) does not yield the same results as defining the directions after taking the closure of the flows, e.g., consider Figure 2. We have seen that x 1 > 0 for all β 1 ∈ R >0 , and thus z + 1 = 1.But in the closure, x 1 = 0 is possible, while still z + 1 = 1 holds, that is, (3) is violated.
• For energy networks β a = ∞ can be interpreted as if the arc is combined with a switch/valve which is turned off/closed.

Acyclic Subgraphs and Computational Complexity
We next obtain a complete description of P AS if the graph is planar by using known results from the literature.Indeed, acyclic (z + , z − ) ∈ {0, 1} 2A correspond to acyclic subgraphs of the digraph D = (V, A).The corresponding acyclic subgraph problem was investigated by Grötschel et al. [16].The acyclic subgraph polytope is the convex hull of incidence vectors of acyclic arc sets in a given digraph.Grötschel et al. showed that for planar graphs a complete description of the acyclic subgraph polytope is given by the variable bounds and so-called dicycle inequalities.These inequalities are based on the set of all dicycles (directed cycles) in D : Note that the anti-parallel arcs {a, a} form a particular dicycle in D .Thus, translated to our setting, we obtain the following: Corollary 10 (Grötschel et al. [16]).If D is planar then Note that planarity is a reasonable assumption for real world physical networks.In general networks, however, optimizing over P AS (in fact, over all four polytopes) is NP-hard.
Lemma 11.Linear optimization over P AS is NP-hard.
Proof.Grötschel et al. [16] already observed that linear optimization over the acyclic digraph polytope is NP-hard, since finding a maximum acyclic subdigraph is NP-hard -this problem is complementary to the feedback arc set, which has been proven to be NP-hard by Karp [21].We note the graph in the reduction is simple and does not contain anti-parallel arcs.When optimizing over P AS and considering D , we can choose the weight to be 0 for either a or a, depending on which direction is present in the original digraph.Thus, optimization over the acyclic subgraph polytope is equivalent to optimization over P AS .
Lemma 12.Given a directed graph with flow bounds and supplies/demands b, it is NPcomplete to decide whether there exists an acyclic b-flow.
Proof.If there exists an acyclic flow, there exists one with polynomial encoding length in the size of the instance.Moreover, acyclicity can be checked in polynomial time.Thus, the problem is in NP.
Consider an instance of the independent set problem: Given an undirected graph G = (V, E) and an integer K, the question is whether there exists an independent subset of nodes of size at least K, i.e., no two selected nodes are connected by an edge.Construct the following directed graph where s and t are two new nodes and v , v are distinct copies of v ∈ V .The arcs in A are constructed as follows: For each edge {u, v} ∈ E we add two arcs (u , v ) and (v , u ); the corresponding flow bounds are such that the flow on these arcs is fixed to 1.Moreover, for each node v ∈ V , we add arcs (v , v ) as well as arcs (s, v ) and (v , t); the flow on these arcs is restricted to lie in Note that because of the flow bounds, the direction of the flows is fixed.However, x a can still be 0 on arcs of type a = (v , v ).
Consider an independent set S ⊆ V of size K in G. Then there exists an acyclic flow in D: For each v ∈ S, the arcs (s, v ), (v , v ), and (v , t) have a flow value of 1.The arcs (s, v ), (v , v ), and (v , t) for v / ∈ S have flow 0. The flow on all other arcs is fixed to 1.It is easy to see that this forms a b-flow.Moreover, it is acyclic.Indeed, because of the flow bounds, the only directed cycles are (u , u , v , v , u ) for an edge {u, v} ∈ E. Since S is independent, the flow on either (v , v ) or (u , u ) is 0. Thus, in each cycle there is at least one arc with zero flow.
Conversely, let x ∈ R A be a feasible acyclic b-flow and define S := {v ∈ V : Because of the demand of −K at t and the flow bounds, we have |S| ≥ K.Moreover, S is independent.Indeed, if there would exist an edge {u, v} ⊆ S, then x would contain a cycle (u , u , v , v , u ).
As a consequence, we cannot expect to obtain tractable linear descriptions for P AS and P AF in general graphs.
Obviously, acyclic subgraphs are not an accurate model for the feasible flow directions, for instance, since proper disconnected subgraphs might not even support a feasible flow.Nevertheless, the acyclic subgraph polytope is well investigated and provides a relaxation.

Acyclic Subgraphs with Sources and Sinks
To obtain a polytope contained in P AS , but closer to P AF , we use the knowledge of sources and sinks in the network.For b ∈ R V , define the set of sources V + := {v ∈ V : b v > 0}, the set of sinks V − := {v ∈ V : b v < 0}, and the inner nodes In the following, for some arc set S ⊆ A we use the shorthand notation z + (S) = a∈S z + a and z − (S) = a∈S z − a .For s ∈ V + , there has to exist at least one arc with flow away from the source s.Similarly, for t ∈ V − , there exists at least one arc with flow towards t.This can be expressed via the valid inequalities Furthermore, for every inner node v ∈ V 0 in-and outflow have to be balanced, because of (2a), that is, there has to exist flow to v if there is flow from v to another node, and vice versa.Thus, if there is an arc (u, v) ∈ A(z + , z − ), there has to exist another node w with (v, w) ∈ A(z + , z − ) as long as z + a = z − a = 0 if x a = 0.There are several possibilities to represent this by linear inequalities.We introduce two options, which differ in their strength and number of added inequalities.
Given a node v ∈ V 0 , the first option is to add an inequality for both directions of every arc incident to v. The inequalities are Here, the first two inequalities imply that if arc a incident to v is oriented away from v, then there has to exist another arc that is oriented towards v.The other two inequalities imply the converse.This discrete representation of flow conservation requires 2 v∈V 0 deg(v) inequalities.
Another option is to aggregate the first two inequalities and the last two inequalities, which yields Again the first inequality implies that if there is an outgoing arc of node v, there has to exist an incoming arc, while the second inequality implies the converse.This representation usually has fewer inequalities: 2 |V 0 | instead of 2 v∈V 0 deg(v).However, while both variants allow for the same integral points, the following example shows that they differ in the strength of their LP-relaxations.
Example 13.Consider again the graph shown in Figure 2. Here, z + 1 = 1, z + 3 = 1 2 , z + 5 = 1, and the remaining variables equal to 0 is feasible for the inequalities (10).In fact, it is a vertex of the LP-relaxation of (7) with the additional constraints ( 8) and (10).However, this solution is not feasible for (9).Furthermore, all feasible solutions of (9) are feasible for (10).This shows that the first option yields tighter LP-relaxations.
Note that for this and subsequent examples we used polymake [1,13] to compute vertices, dimension, affine hull, and facets of polytopes.
In the following we will therefore concentrate on (9) and investigate the formulation given by these inequalities as well as the dicycle inequalities We define the polytope of acyclic flows with sources and sinks as ) is feasible for ( 8), ( 9), (11) .
The model derived here turns out to be equivalent to the one investigated by Becker and Hiller [19,2,3], which can be seen using the analysis in the next section.

Analysis of Acyclic Subgraphs with Sources and Sinks
In this section we analyze the polytope P AS± .We start by proving a key insight, which helps to derive several results in the following.Note that we define paths to be simple, i.e., no node appears twice in the path.We call a directed path a source-sink-path if it starts in V + and ends in V − .Proposition 14.Let (z + , z − ) ∈ P AS± ∩ {0, 1} 2A .For every arc in A(z + , z − ), there exists a directed source-sink-path containing this arc.In particular, D(z + , z − ) contains at least one path leaving each node in V + and at least one path entering each node of V − .Moreover, if b = 0, then P AS± = {0}.
Proof.Consider an arbitrary arc a = (u, v) ∈ A(z + , z − ).We construct a path from v to V − and a path from V + to u, which together with a yield the desired path.Note that these three paths are necessarily simple, since otherwise A(z + , z − ) contains a cycle.
If b = 0, by (8a) there exists at least one path leaving each node in V + and by (8b) there exists at least one path entering each node in V − .
Finally, let b = 0. Then there are no sources and sinks, i.e., the construction above would either terminate at a node with degree 1 or produce a cycle, which contradicts either (9) or (11).Thus, We obtain the following first consequence: Then, if there are no flow bounds, Proof.In the case b = 0, we have P AS± = P AF = {0}.Hence, we only have to consider the case b = 0.
We first show that P AF ⊆ P AS± .(Note that we do not need any assumption on the flow bounds or on the sizes of V + and V − for this direction and that it suffices to prove the inclusion for all integer points.)Let (z + , z − ) ∈ P AF ∩ {0, 1} 2A with corresponding acyclic b-flow x ∈ R A .Since b = 0, x is nonzero.This implies that there is at least one path with nonzero flow leaving each node in V + and at least one path with nonzero flow entering each node in V − .Thus, constraints (8) are satisfied.Due to flow conservation, (9a) -(9d) hold.Further, since x is acyclic (11) holds.
We continue by proving the other direction, i.e., P AS± ⊆ P AF .The proof works for the single sink case V − = {t}, since the single source case is analogous.Let V + = {s 1 , . . ., s k }, ) and (t 1 , t 2 ) has to be positive.Nevertheless, if b t 1 = 0, they need not be used in P AS± , i.e., P AS± ⊆ P AF .and let (z + , z − ) ∈ P AS± be a binary vector.We first construct an acyclic flow , and x a = 0 otherwise.By scaling, we will obtain a b-flow.Start with the zero-flow x = 0 and define P 1 = • • • = P k = ∅.We then pick an arc a ∈ A with either z + a = 1 or z − a = 1 and x a = 0. Then by Proposition 14 there exists a path P from a source s i to the sink t in D(z + , z − ) containing a (if z + a = 1) or a (if z − a = 1).Add P to the set P i and augment x along P by one unit, by increasing x a by 1 for every arc a ∈ A ∩ P and decreasing x a by 1 for every arc a ∈ A such that a ∈ P .For a ∈ A define ∆(P ) a = 1 if a ∈ P and ∆(P ) a = −1 if a ∈ P and ∆(P ) a = 0 otherwise.Then the new flow is x + ∆(P ).Since A(z + , z − ) does not contain both a and a , x a can only be increased (if a ∈ A(z + , z − )) or decreased (if a ∈ A(z + , z − )) by augmenting flow along another path, even for another source-sink pair.Therefore, we can iterate augmenting flow for the remaining arcs with no flow and thereby construct a flow with the desired flow directions.Note that We still have to scale the flow such that it is a b-flow.First consider each i ∈ {1, . . ., k} in which P i = ∅.By Proposition 14 there exists an s i -t-path P in D(z + , z − ).Then set P i = {P }.We now define the scaled flow Since b s i and |P i | > 0, the scaling is valid, does not change flow directions, and yields a b-flow.Remark 16. Figure 3 shows that Theorem 15 does not hold with more than one source and sink.
Proposition 14 also helps to determine the structure of integer points in P AS± .For a subset S ⊆ A, let χ(S) be the incidence vector Corollary 17.Each integer point (z + , z − ) ∈ P AS± is the incidence vector of a union of source-sink-paths in D that does not contain cycles and conversely.
Proof.For each a ∈ A(z + , z − ) there exists a source-sink-path P a in A(z + , z − ) that contains a by Proposition 14.Then each arc in A(z + , z − ) is covered by ∪ a P a and the union does not contain cycles.Thus (z + , z − ) = χ(∪ a P a ).
Conversely, the incidence vectors of a union of source-sink-paths that does not contain cycles clearly satisfies ( 8), (9), as well as (11) and is therefore contained in P AS± .
Note that the union of source-sink-paths can contain cycles, even in the single source and sink case, see the example in Figure 1.
Another consequence of Proposition 14 is that z + a and z − a can be fixed to 0 or 1 in some cases.We first need the following definition.Recall that we assume that D is weakly connected and consider two arcs a 1 and a 2 ∈ D such that neither is a bridge and D − {a 1 , a 2 } has exactly two weakly connected components D 1 and D 2 .Then {a 1 , a 2 } is called a cut-pair.Assume that D 2 contains neither source nor sink and that a 1 enters and a 2 leaves D 2 (by reorientation).We call D 2 input-output subgraph.Note that D 2 might have no arcs, in which case a 1 and a 2 form a directed path.
Lemma 18.Let D be the given weakly connected, simple digraph with sources V + and sinks V − .Then the following holds for every (z + , z − ) ∈ P AS± .
1.If there is no source-sink-path in D containing a ∈ A ( a ∈ A) then z + a = 0 (z − a = 0).2. Let a be a bridge, i.e., D − a is not weakly connected.Then the following holds for each of the two connected components induced by B ⊆ V.
Let there exist an input-output subgraph of D with entering arc a and leaving arc a .Then z + a = z + a and z − a = z − a .Proof.In all cases, it suffices to consider integer points (z + , z − ) ∈ P AS± , since the statement then holds for the convex hull P AS± .1. Suppose that there is no source-sink-path in D containing a ( a).By Proposition 14, A(z + , z − ) cannot contain a ( a).Thus, z + a = 0 (z − a = 0).2. In case (2a), B contains neither a source nor a sink.Since the paths are simple and D(z + , z − ) is acyclic, no path can enter B. In case (2b), at least one source-sink-path has to enter B. In case (2c), at least one source-sink-path has to leave B.
3. Every union of source-sink-paths using a 1 also has to use a 2 .This implies the given equations.
1.The results in Lemma 18 are similar to the ones by Becker and Hiller [2,3], but in a different notation.
2. The existence of a node of degree 2 is a special case of Part 3 of Lemma 18.

3.
It is an open question whether the conditions of Lemma 24 define the affine hull of P AS± ; but see Proposition 28.
A natural question is how the conditions in Lemma 24 can be checked.It turns out that the condition of Part 1 is hard to check, even for the single source and sink case.
Proposition 20.Given a directed graph with source node s, sink node t and some arc a = (u, v), it is NP-complete to decide whether there exists a (simple) s-t-path that contains a.
Proof.Consider the k-vertex disjoint paths problem, which consists of finding vertex disjoint paths from s i to t i for a given set of node pairs (s 1 , t 1 ), . . ., (s k , t k ).Obviously, finding an s-t-path that uses a is the special case of finding 2-vertex disjoint paths between (s, u) and (v, t).Fortune et.al [10] proved that the vertex disjoint paths problem on general directed graphs is NP-complete even for fixed k ≥ 2. Corollary 21.Linear optimization over P AS± is NP-hard, even if there is a single source s and sink t.
Proof.Consider the linear function that maximizes z + a for some arc a over P AS± .The optimal value is 1 if and only if there exists an s-t-path through a.The results then follows by NP-hardness of determining the latter by Proposition 20.
Remark 22. Schrijver [35] showed that for fixed k and planar graphs, an s-t-path that contains a given arc can be found in polynomial time.Recently, Fakcharoenphol et al. [9] showed that one can compute the set of all arcs that are not contained in an s-t-path of a planar graph in linear time.
Remark 23.With respect to Parts 2 and 3 of Lemma 18 the following holds.Bridges in graphs can be found in linear time, see Tarjan [39].Moreover, checking whether a source and sink are in the same connected component can be done by breadth-first search in linear time, see, e.g., Korte and Vygen [23].Moreover, after bridges have been removed, the linear time algorithm of Mehlhorn et al. [28] outputs a cut-pair if one exists.More input-output subgraphs can be produced using this algorithm iteratively.

Analysis of the Single-Source and Sink Case
In this section, we provide a further analysis for the special case of a single source s and sink t.This implies that the balanced b ∈ R V satisfies b s = −b t ≥ 0 and b v = 0 for all v ∈ V \ {s, t}.To simplify notation, we orient the arcs incident to the source s and sink t such that δ − (s) = ∅ and δ + (t) = ∅ holds.
Lemma 24.Let D be the given weakly connected, simple digraph with source s and sink t.Then for every (z + , z − ) ∈ P AS± , z − a = 0 holds for every a ∈ δ + (s) and for every a ∈ δ − (t).
Proof.First consider an integer point (z + , z − ) ∈ P AS± .Assume the statement does not hold, and let a = (s, v) with z − a = 1.Then going backwards from v similarly to the proof of Proposition 14, there exists an s-v-path.This would close a cycle, hence, z − a = 0 holds for all a ∈ δ + (s).For a ∈ δ − (t) we can argue analogously.Since P AS± is the convex hull of the integer points, the statement holds.
Remark 25.Let S ⊂ V with s ∈ S, t / ∈ S and consider the s-t-cut inequalities Because of Corollary 17 these inequalities are valid for all integer points in P AS± and thus for their convex hull.However, they are weaker than the inequalities (8) and (9).Indeed, the s-t-cut inequalities together with nonnegativity provide a complete linear description of the dominant of the s-t-path polytope, see, e.g., Schrijver [36,Thm. 13.1].Moreover, as an example consider the incidence vector of the union of (at least) one s-t-path and some node-disjoint disjoint arc.This vector is feasible for the s-t-cut inequalities, but not for (8) and (9).However, these inequalities can be strengthened as follows.
Lemma 26.Let D be a simple connected digraph with source s and sink t.Let S ⊂ V with s ∈ S, t / ∈ S. Then the inequalities are valid for P AS± .
Proof.Inequality (13a) is satisfied by all solutions in P AS± with z − a = 0 by Remark 25.Let (z + , z − ) ∈ P AS± be an integral solution with z − a = 1.By Proposition 14, D(z + , z − ) contains an s-t-path P that uses a. Then P has to cross the cut at least twice and therefore z + (δ + (S)) + z − (δ − (S)) ≥ 2. Thus, (13a) holds for the convex hull of these integer points, i.e., P AS± .The validity of (13b) can be seen similarly.
We can prove the following concerning the affine hull of P AS± .Proposition 28.Assume that in D there exist two arc-disjoint s-t-paths and for every arc a ∈ A \ { a : a ∈ δ + (s) ∪ δ − (t)}, there exists two s-t-paths that are arc-disjoint except for a.
where δ + , δ − and deg are with respect to the original digraph D.
Consider an arbitrary arc a ∈ A \ (δ − (s) ∪ δ + (t)).By assumption there exist two s-tpaths P 1 and P 2 that are arc-disjoint, except for a.Let P1 := P \ {a} and P2 := P \ {a}.Then for i ∈ {1, 2}: Adding the first two equations yields This implies c a = 0, and therefore no equation can be active for a.Thus, the only possible equations are the ones of Lemma 24.Together, they reduce the dimension by |δ − (s)| + |δ + (t)|, which shows the claim.
Note that the assumptions of Proposition 28 rule out bridges as well as input-output subgraphs.Moreover, note that the assumptions also rule out the existence of an s-t-arc, since the only s-t-path using this arc would be the arc itself.Nevertheless, we could relax this condition for an s-t-arc, but then, if (s, t) ∈ A, the formula would have to be dim P AS± = | A| − deg(s) − deg(t) + 1 due to double counting.

Model Implementation
Recall from Example 1 that in stationary gas transport, the potentials are the squares of the pressures at the nodes.Since some of the models of network elements other than pipes cannot be (linearly) expressed in terms of the potentials, our model includes pressure variables p u for all nodes u ∈ V, which are coupled with the potentials through the equation p 2 u = π u .To avoid unnecessary nonlinear equations, our model only contains the potential variables, where they are actually needed.Note that this was a modeling choice and one could also only introduce the pressure variables as needed.
• Pipes are handled in the way described in Example 1 with ψ(x a ) = |x a | x a and the resistances given by , where L a and D a are the length and diameter of the pipe a = (u, v) ∈ A, respectively, R s is the specific gas constant, T m is the temperature (assumed to be constant), z m is the constant z-factor, and λ a the friction coefficient using the formula of Nikuradse [29,30], where k a is the roughness of the pipe.The z-factor is computed using the formula of the American Gas Association see Králik et al. [24], where p c and T c denote the pseudocritical pressure and pseudocritical temperature of the gas.Then z m := z(p m , T m ) with the mean pressure value derived from the pressure bounds at nodes u and v.
• Valves control the gas flow.For open valves, the pressures on both sides are equal.For closed valves, the flow is 0 and the pressures are decoupled.These conditions are modeled using a binary variable in a straight-forward way.
• Compressors can increase gas pressure.Several compressors are often connected by piping and valves in compressor stations.We approximate the operation states of such stations by polyhedra in terms of flow, input and output pressure, see Hiller and Walther [20].The operation modes of compressors are modeled using binary variables: the compressor can be turned on/off or it can be in bypass, i.e., flow can bypass the compressor or can flow in the opposite direction.The possible modes are: compressor on/bypass closed, compressor off/bypass open and compressor off/bypass closed.
• As objective function we chose the maximization of the sum of pressures.Several other objective choices exist, see, e.g., [17].
For handling acyclic flows we add the variables z + , z − which are coupled with the flows by using the linear relaxation (5) instead of (3).For all pipes a = (u, v) ∈ A they are coupled with the pressure variables by inequality Moreover, they can be integrated in some of the other network elements in a straightforward way, e.g., binary variables for valves are coupled with the direction variables, such that z + a = z − a = 0 if the valve is closed.Depending on the model variant, we add the binary representation of flow conservation given by the inequalities ( 8) and (9).The dicycle inequalities (11) can be added for a cycle basis or for all cycles.Therefore, a cycle basis is computed as follows: after computing a spanning tree by breadth-first-search, each non-tree arc defines a cycle in our cycle basis.The possible combinations of these cycles are then enumerated.Two cycles can be combined into a new cycle, if their symmetric difference induces a connected subgraph, where all nodes have degree two.Note that the enumeration and checking if two cycles can be combined, takes only a fraction of a second for the networks considered here.
Note that we forbid flow in cycles over compressor stations, because such cycles would significantly increase the temperature of the gas.This could only be controlled in transient models where the temperature of the gas is kept under control as well.

Results
In order to test the effects of the different conditions we performed computations with the following model variants: FLC+AC as variant FLC plus dicycle inequalities for all cycles.
The computations were performed on a cluster with 3.5 GHz Intel Xeon E5-1620 Quad-Core CPUs, having 32 GB main memory and 10 MB cache running Linux.We used SCIP-7.0.0 [11,37] with a one hour time limit and we used CPLEX-12.10.0 as LP-solver.
Since the GasLib-40 instance is rather small, the solving time for all models is less than a second (and hence not reported in detail).Thus, this does not allow to draw conclusions on the different model variants.Nevertheless, the instance gives insight on some advantages of using the flow direction variables.Table 1 shows statistics on the flow bounds of the pipes after presolving.Column "#fixed flows" shows the number of pipes with fixed flow, column "#fixed dirs" the number of pipes with fixed flow direction, and columns "x mean " and "x mean " the mean arithmetic lower and upper bounds of the flows.While the number of fixed flows is the same for all models, and mainly depends on the graph and the position of sources and sinks, the number of fixed directions and the mean arithmetic flow bounds can be improved by using the flow direction variables.This is also illustrated by Figure 6, which compares   [14] and the references therein.Indeed, Becker and Hiller used OBBT [3] to further strengthen flow bounds.
The above results show that the mean arithmetic flow interval for model NFD is more than four times as large as the flow interval of model FLC+AC.Since tighter variable bounds typically lead to smaller branch-and-bound trees, this positive effect on the flow bounds will also be reflected in the solving times of larges instances as can be seen by the following results for the GasLib-582 network.
The results of each model variant on the GasLib-582 network are given in Tables 2, 3, and 4 aggregated over all 4227 scenarios.In Table 2, column "optimal" gives the number of feasible scenarios solved to optimality, "feasible" the number of scenarios for which a feasible solution could be found, but could not be solved to optimality, "limit" the number of scenarios running into the time limit, "infeasible" the total number of scenarios that have been determined to be infeasible, and "inf-presol" the number of instances that have been determined to be infeasible during presolving.
Table 3 provides statistics on the (geometric mean) running times of the model variants.Here, column "to optimality" is the mean geometric time it took to prove optimality, the column "to first" gives the mean geometric time until the first feasible solution was found, "infeasible" the geometric time to proof infeasibility, "totalgeometric" the total geometric time, and "totaltime" show the total computational time in hours.
In Table 4 similar statistics to Table 1 are given.Here, the results are averaged over all scenarios.Additionally, the columns "#min dirs" and "#max dirs" show the minimal and maximal number of fixed flow directions over all scenarios.Remark 32.The numerical results shown in Table 2 are consistent in the sense that all variants identify the same 2179 infeasible scenarios.Moreover, for all other scenarios feasible solutions were found by at least one model variant.In fact, only 20 scenarios could not be solved to optimality.For these instances, at least one model variant found a feasible solution with optimality gap 0.4 % or better.
The results clearly show that determining infeasibility seems to be easy in most cases.With all model variants, almost all infeasible scenarios could already be identified during presolving.That is, our acyclic flow models only slightly improve the computations here.The numbers for feasible instances, however, show a completely different picture.Comparing the solving times in Table 3 for the basic model without any binary direction variables (i.e., NFD) with the model enhanced by P AS± (i.e., FLC+AC) shows a speed-up factor of ∼ 4.9 in geometric time to prove optimality, and a speed-up factor of ∼ 7.2 in the total running times.Moreover, Table 2 shows that with model NFD for almost half of the feasible instances no feasible solution has been found, while with model FLC+AC almost all feasible instances could be solved to optimality.
A partial explanation for the performance improvement is as follows.With the flow direction variables and the additional constraints we represent properties of feasible solutions, which are otherwise not included in the initial model/relaxation.Moreover, the heuristics and presolving techniques can detect more variable fixings, implications and reductions based on the flow conservation and the flow direction variables.For example, in diving heuristics, it is easier to detect infeasibility based on the binary variables without having to consider the nonlinear physics.This then leads to tighter variable bounds (after presolving), which can be seen for the flow variables in Tables 1 and 4. Since having tight variables bounds is important to derive good relaxations for the nonlinearities, the LP-relaxations are already stronger early in the branch-and-bound tree.Moreover, tighter variable bounds (typically) lead to smaller branch-and-bound trees, since the search space is smaller.Indeed, we can observe this effect for the computations on the network GasLib582.The arithmetic mean number of branch-and-bound nodes (rounded up) for the optimally solved (feasible) scenarios with model variant NFD is 168 970, while it is 35 330 with model FLC+AC.That is, in model NFD it takes about 4.8-times as many nodes in comparison with variant FLC+AC, which is almost the same ratio as the speed-up of the geometric solving times for these scenarios.[17].Moreover, we also used the model FLC+AC to derive the computational results there.In particular, the results for the ODE model show an even stronger effect of using FLC+AC; see Tables 5 and 6.Note that for these computations, OBBT was used to strengthen bounds of the flow variables.With model NFD, almost half of the scenarios ran into the time limit without a feasible solution.Only 41 scenarios could be solved optimally and feasible solutions were found for only 93 further scenarios.In contrast, with model FLC+AC all scenarios were either proven to be infeasible or a feasible solution was found.Moreover, 84 % of the feasible instances were solved optimally.

Conclusions and Open Questions
In this article we have investigated the usage of the property that potential-based flows are acyclic.We derived combinatorial models and investigated their properties.One of the main goals was to use these models to speed-up optimization problems over potential-based flows.Indeed, our computational results show a speed-up of at least a factor of 3 if cycle information and knowledge on sources and sinks is used.The time to prove optimality shows a speed-up factor of almost 5, and the total running time speed-up of about 7.
In our computations we added the dicycle inequalities up front.It would be interesting to compare the performance of an algorithm that dynamically separates dicycle inequalities.
There are several open questions that we plan to investigate in the future.One open question is whether the conditions in Lemma 24 suffice to define the dimension of P AS± .Obtaining facets of P AS± would be interesting, in order to possibly further strengthen the formulation of potential-based flows.Such facets might be transferred from the acyclic subgraph polytope.This seems to be difficult, though.Moreover, it would be interesting to investigate whether the property of unique flows could be exploited in combinatorial models as well.

Figure 1 :
Figure 1: An s-t-flow which can be decomposed into two paths (indicated by dashed/dotted arcs) each with flow value 1, but is not acyclic.

Figure 2 :
Figure 2: A diamond shaped network, with source s and sink t and arc labels 1 to 5.

Figure 3 :
Figure 3: The graph shows that Theorem 15 does not hold for |V + |, |V − | ≥ 2: Suppose that b s 2 , b t 2 = 0 and b s 2 < −b t 2 , then the flow on at least one of the arcs (s 1 , s 2) and (t 1 , t 2 ) has to be positive.Nevertheless, if b t 1 = 0, they need not be used in P AS± , i.e., P AS± ⊆ P AF .

Figure 4 :
Figure 4: The graph D(z + , z − ) associated with a vertex of P AS± LP applied to the network of Figure 2.

NFD
no binary variables to represent flow directions; AC dicycle inequalities (11) for all cycles; FLC binary flow conservation (8) and (9); FLC+CB as variant FLC plus dicycle inequalities for a cycle basis;

Figure 6 :
Figure 6: The presolved network GasLib-40 corresponding to variants NFD (left) and FLC+AC (right).The scenario has 3 sources (diamonds) and 29 sinks (circles).Pipes with fixed flow are depicted by , fixed flow directions are shown by , and the remaining pipes (with unfixed flows/directions) are dashed.

Table 1 :
Statistics on the flow bounds after presolving of 39 pipes in network GasLib-40 for all model variants.

Table 2 :
Aggregated results for GasLib582 scenarios for all model variants.

Table 3 :
Geometric means of solving times in seconds and total run time in hours for all model variants for the GasLib582 scenarios.

Table 4 :
Statistics on the flow bounds after presolving of 278 pipes in network GasLib-582 for all model variants.

Table 5 :
Aggregated results for GasLib582 scenarios for the ODE model.

Table 6 :
Geometric means of solving times in seconds and total run time in hours for the GasLib582 scenarios for the ODE model.Remark 33.The fact that potential-based flows are acyclic does not only hold for this algebraic model for stationary gas flows.It also applies to the stationary model based on ordinary differential equations used in