Maximum Throughput Problem in Dissipative Flow Networks with Application to Natural Gas Systems
Abstract
We consider a dissipative flow network that obeys the standard linear nodal flow conservation, and where flows on edges are driven by potential difference between adjacent nodes. We show that in the case when the flow is a monotonically increasing function of the potential difference, solution of the network flow equations is unique and can be equivalently recast as the solution of a strictly convex optimization problem. We also analyze the maximum throughput problem on such networks seeking to maximize the amount of flow that can be delivered to the loads while satisfying bounds on the node potentials. When the dissipation function is differentiable we develop a representation of the maximum throughput problem in the form of a twice differentiable biconvex optimization problem exploiting the variational representation of the network flow equations. In the process we prove a special case of a certain monotonicity property of dissipative flow networks. When the dissipation function follows a power law with exponent greater than one, we suggest a mixed integer convex relaxation of the maximum throughput problem. Finally, we illustrate application of these general results to balanced, i.e. steady, natural gas networks also validating the theory results through simulations on a test case.
I Introduction
We consider a network on an undirected graph , where is the set of nodes with and is the set of oriented edges with . For a dissipative flow network, we have two sets of variables describing the network. Each is assigned a nominal direction and is characterized by a flow from node to node , and each node is characterized by a node potential . Such networks which have a description via potentials are determined by a set of two equations:“flow conservation” and “potential drop relation.” Perhaps the most famous example of the dissipative flow networks is a linear DC resistive electric circuit, where the flows are currents and node potentials are voltage potentials. In this case, the flow conservation equations correspond to the socalled Kirchoff’s first law, and the potential drop relation is the Kirchoff’s second law, or Ohm’s law. Resistive circuits are linear dissipative networks, where the potential drop relation given by the Ohm’s law is linear. In this paper we consider more general dissipative networks where the potential drop relations is a general nonlinear relation described by a monotonically increasing dissipation function.
We formulate a variational representation of the network flow equations in the form of a strictly convex optimization problem, thus proving the uniqueness of the solution. Additionally, this solution can be obtained efficiently by using convex optimization algorithms such as the interior point method. The energy function method for gas flow networks was studied in [1, 2]. The maximum throughput problem is an optimization problem that aims at maximizing the amount of flow delivered from sources to loads, while satisfying bounds on the potentials and possibly other control variables. For general dissipative flow networks, this problem is nonconvex. In the case when the dissipation function is differentiable, we develop a twice differentiable optimization formulation of the maximum throughput problem that is biconvex in the optimization variables in a certain regime. We obtain several properties of the primal and dual form of the energy function that are useful to provide gradient and hessian information to optimization software during numerical implementation. In the process we prove a special case of the socalled “monotonicity property” of dissipative flow networks, also discussed in a companion paper [3].
When the dissipation function takes the form of a power law with exponent greater than one, we provide a mixed integer convex relaxation of the maximum throughput problem, where the degree of nonlinearity is governed by the exponent. Together, we have a framework to obtain a feasible solution using the energy function formulation, and then check the quality of the solution (degree of suboptimality) using the mixed integer convex formulation. As an example, we demonstrate how the setting applies to the natural gas networks and also illustrate its performance on a numerical example.
The rest of the manuscript is organized as follows. In Section II we define dissipative flow networks and equations that govern them. Section III states the maximum throughput problem for a dissipative flow network. Section VIA introduces variational formulations of the network flow equations and their properties and provides a reformulation of the maximum throughput problem using the energy functions. In Section V we introduce a mixed integer convex relaxation for the maximum throughput problem where the dissipation function follows a power law with exponent greater than one. In Section VI we briefly describe natural gas networks and adopt the techniques of the previous Sections to this special enabling case. Finally, in Section VII we provide computational results on natural gas network models and then summarise and present path forward in Section VIII.
Ii Dissipative Flow Networks
A dissipative flow network, where flows on edges are driven by potential difference between adjacent nodes, is described by the following two sets of the Network Flow (NF) equations:

Flow conservation:
(1) where is the injection at node .

Potential drop relation:
(2) The function is a differentiable, strictly monotonically increasing odd function which is called the dissipation function.
Several networks can be described within the framework. For example, DC resistive electric circuits [4], AC power flows in the lossless approximation [5, 6, 7, 8], balanced natural gas networks [9, 10, 2, 11], traffic flows [12, 13, 14] etc. The term allows us to model additional network elements in the actual networks such as compressors in gas networks, voltage transformers in DC resistive circuits and phase shifters in the approximate AC power flows.
The system of potential drop relations (2) is translation invariant, i.e., if is a solution of the NF Eqs. (1,2), then a translation of the potentials by a constant, , is also a solution. To remove this degree of freedom, the potential at a particular node (node ) is fixed at a certain value . In the case of the resistive circuits, this node is called “ground”, in AC power system respective node is called “slack bus”, and this term is also adopted to the case of the gas networks, where the slack bus (where the pressure is fixed) is usually assigned to the largest source of gas in the system. Assuming that parameters are known, we are now left with the number of variables, (), equal to the number of equations. However, this still does not guarantee existence and uniqueness of the NF Eqs. (1,2). But as we will see shortly, in Section IIA, this is true for dissipative networks with a strictly monotonically increasing dissipation function. Further, observe that the set of flow conservation equations (1) automatically implies total flow balance
(3) 
Iia Primal Energy Function and Uniqueness of the Network Flow Solution
Let be the antiderivative of , i.e, . The following useful remark follows immediately from the fact that is monotonically increasing.
Remark 1
The function is convex.
The set of NF equations (1) and (2) can be rewritten in the variational form as the optimal solution of the following convex optimization problem [1, 2]:
(4)  
s.t.  (5) 
The optimization problem in (4)(5) is clearly convex, since by Remark 1 the cost function is convex, and the constraints are linear. To see why the solution to (4)(5) yields the NF equations, we can first form the Lagrangian as
(6) 
where are the dual variables corresponding to the constraints (5). To minimize the Lagrangian, we take its derivative w.r.t. to zero to get
(7) 
which means that the dual variables serve the role of the node potentials. Since the optimization problem in (4)(5) is strictly convex, we also conclude that the NF Eqs. (1) and (2) have a unique solution.
IiB Dual Energy Function Representation
We now present another variational formulation of the NF equations via the dual of the optimization problem in (4)(5). Since is monotonically increasing, it is invertible. Let denote the inverse of . Let be the antiderivative of . The dual of (4)(5) can be expressed in terms of the node potential variables as:
(8) 
The above problem, when stated as a minimization, is a convex optimization problem. Since the Slater’s condition [15] always holds for an underdetermined system of linear equality constraints, we have strong duality and the optimal value of (4)(5) and (8) are the same.
Iii The Maximum Throughput Problem
In a dissipative network the nodes with are called sources and those with are called sinks. Assume that there are constraints on the minimum and maximum value of the allowed node potentials, , in the form
(9) 
The maximum throughput problem seeks to maximize the amount of flow that is delivered to the loads while still satisfying the box constraints on node potentials in (9). Formally:
Maximum Throughput in Dissipative Flow Networks:
(10)  
s.t.  (11)  
(12)  
(13) 
The potential drop relations (12) are in general nonconvex and make the maximum throughput problem a nonconvex problem.
Iv Representing the maximum throughput problem using energy functions
In this Section, we provide another alternative (but, obviously, equivalent) formulation for the maximum throughput problem using the variational representation of the NF equations in terms of the energy function. To achieve this goal we, first, discuss optimal solution and value at the optima of the energy function.
Iva Properties of the energy functions
Define the following convex function given by
(14) 
and let
(15) 
The function is the convex conjugate of with respect to . Analogously, define
(16) 
and
(17)  
s.t.  (18) 
The following remark follows directly from the definitions above.
Remark 2
Lemma 1
The following statements follow from standard results of the convex conjugate and the duality theory [16].

The Hessian of can be computed with respect to as follows
(21) where
(22) (23) (24) (25) 
The inequality
(26) has exactly one solution given by .

The statement in Eq. (19) follows directly from the theory of convex conjugates. The proof of (20) is similar, but we include it here for completeness. Recall that . So, it is enough to prove that . By definition of , we have
So,
Using Eq. (7), one derives
(27) (28) Since , we have that . The proof is complete by substituting this result in Eq. (28).

Follows from standard properties of convex conjugates and direct computations.

By definition is the minimum of w.r.t. . Since is strictly convex, only its unique minima can satisfy (26).
For the statement in Lemma 1 (b) to be welldefined, we need that the matrix is invertible. Note that since the pressure at node is fixed at , indices in Eq. (23) and Eq. (25) run from to . From the expressions for the elements of the hessian, we derive that is a submatrix of a weighted graph Laplacian matrix with edge weights given by . Since the graph is assumed to be connected, it follows that has full rank and hence is invertible [17], [18].
Let us now take a slight detour from the main goal (reformulation of the maximum throughput problem) and state a certain “monotonicity property” that holds for dissipative FN. For a more general statement and proof of the monotonicity property without assuming the differentiability of the dissipation function as well as its implication on robust optimization on dissipative flow networks, we refer to [3].
Corollary 1
(Monotonicity Property) Consider a dissipative flow network where the parameters and the potential are fixed. Let and be two sets of injections. The injection at node is automatically determined by the total flow balance Eq. (3). Assume that for all one has . Let and be the corresponding potentials resulting from the two injection profiles. Then one arrives at for all .
IvB Equivalent formulation of the maximum throughput problem
Equipped with Lemma 1, one restates the maximum throughput problem (10)(13) as
(29)  
s.t.  (30)  
(31) 
This formulation is exactly equivalent to the original maximum throughput formulation. Let us present two other alternative formulations, which (as we see below) may have some practical advantages.
Energy function formulation :
(32)  
s.t.  
Energy function formulation :
(33)  
s.t. 
In (32) is a small constant. In (33) is a large positive constant that helps to enforce the constraint (30).
Eq. (32) and Eq. (33) are useful because, first, the two expressions are twice differentiable and hence allow application of second order gradient descent methods for local optimization. Second, when is not an optimization variable, the expressions are biconvex in and which means that we can use here standard heuristics, such as alternating minimization over and . Moreover, solvers using second order methods can be provided with first and second derivative information by using Eqs. (19)(20) and Eqs. (23)(25). The biconvexity is lost when is also subject to optimization, but the twice differentiability remains.
V Mixed Integer Convex Relaxation
In this Section we formulate a mixed integer convex relaxation of the maximum throughput problem (10)(13) in the special case when the dissipation function is a power law, i.e., with . Let us start introducing binary variables associated with each edge that denote the flow direction, i.e., . We can then restate the nonconvex constraint (12) as
(34)  
or equivaliently  (35) 
Next, we apply the McCormick relaxation [19] to the bilinear term on the lhs and relax the equality to an inequality to get the following pair of inequality relaxations of (12):
(36)  
(37) 
Since the constraints (36) and (37) are mixed integer convex constraints. The McCormick relaxation is tight when at least one of the variables is binary, the only relaxation in this set of constraints consists in replacing the equality in (34) by an inequality.
The mixed integer convex relaxation along with the formulation in Section VIA provide the framework for finding a feasible solution of the maximum throughput problem as well as for testing the quality of the solution through comparison with the lower bound.
Vi Natural Gas Networks
In this Section, we introduce natural gas networks and their properties and apply the techniques of Section VIA and V to the maximum throughput problem over these networks. Natural gas has emerged as an alternative to other fuel sources such as coal and fuel oil due to its lower carbon footprint. Improved technology for tapping new natural gas reserves has increased the availability and reduced its cost leading to increasing utilization of natural gas in power generation (gas turbines) and household heating. Natural gas networks use a system of high pressure transmission pipelines to transport gas from sources to loads. The flow of gas in pipelines is driven by gradient in pressure along the length of the pipe. Pressure falls rapidly with distance along the pipelines, so compressor stations are installed to locally boost the pressure and maintain the desired pressure and throughput at the load end. Typically, compressor stations are placed every 50100 miles. Many of these compressors are powered by a fraction of the gas in the transmission line itself. The fraction of gas consumed can vary from . Due to the increasing demand for natural gas, there is an increasing need for optimization of pipeline operational planning. In what follows, we describe two such optimization problems of interest.
The first is the problem of minimizing the cost associated with compression. In the literature, study of this problem can be traced back to [20] (see also [21] for a review of this problem) where a Dynamic Programming approach was used to solve it. More recently in [22], a convex Geometric Programming (GP) based approach was developed for compressor cost optimization in tree networks. This approach offers several desirable properties such as fast convergence, and scope for generalization to a heuristic for loopy networks. Since the compressors consume a part of the gas itself, the compressor cost optimization problem is also related to the next problem below.
The second is the maximum throughput problem introduced in Section III. This is particularly relevant in, for example, Norway, which produces more gas than can be domestically consumed. This situation creates strong economic incentives to sell the excess gas to the rest of Europe. In this case the primary objective of the gas pipeline operator becomes to maximize the throughput of the network within the engineering and contractual constraints that must be satisfied.
In the rest of this Section, we describe the equations governing the flow of natural gas over natural gas networks in the balanced (steady state) regime and apply the machinery of Section VIA and V to this special case.
Natural gas networks in the steady state are governed by the following set of NF equations:
(38)  
(39) 
where is the flux of gas (flow per unit length) from node to node in the steadystate, and is the pressure at node . The quantity models the pressure boost provided by the compressor on the edge . For a more detailed derivation of the steadystate gas flow equations from dynamic fluid mechanic equations, we refer to [22].
From Eqs. (38,39) it is apparent that natural gas networks in steady state can be represented within the framework of dissipative flow networks introduced in Section II by setting
(40)  
(41) 
Via Energy function representation for maximum throughput in natural gas networks
In this Section, we adopt the reformulation of the maximum throughput problem from Section to the case of the natural gas networks. The expressions for functions and in Eqs. (14,16) for gas networks are given by
(42)  
(43) 
The inverse is given by
(44) 
and the hessian matrix is given by
(45)  
(46) 
The energy function formulation of the maximum throughput problem for the gas networks can be obtained substituting the expressions above to Eqs. (2931).
ViB Mixed Integer Quadratic Programming (MIQP) relaxation of the maximum throughput in the natural gas networks
The gas network dissipation function (39) follows a power law with . So the Mixed Integer Convex relaxation of Section V is now a Mixed Integer Quadratic Program (MIQP). We write the full program below for completeness.
MIQP for Maximum Throughput in natural gas networks
(47)  
s.t.  (48)  
(49)  
(50)  
(51)  
(52)  
(53) 
As remarked before, the only relaxation in the above program is the conversion of equality into inequality in Eqs. (49,50). This corresponds to allowing “decompression”, i.e., the pressure drop is now allowed to be larger than that implied by the flow. A similar relaxation was used in the Geometric Programming approach in [22] to preserve convexity, and was observed to perform well in numerical experiments.
Vii Numerical Computations
In this Section, we test the optimization formulations of Sections VIA, V on a model of natural gas network. Fig. 1 illustrates this exemplary network. There are nodes and edges. The nodes marked with green are sources and those with red are sinks. The arrows represent edges and the direction of the arrow denotes the nominal orientation of the edge. The cost is set to zero for sources and nonconsumers. For consumers, we set proportional to the nominal demand at node . We consider three different sets of upper bounds for the pressure squared variables : . In all cases we set . The slack node is set to the maximum pressure.
We solve the maximum throughput problem using the energy function formulation (29)(31) and the MIQP formulation (47)(53) both with and without compressors included into the set of optimization variables. To test the quality of the solution and bounds thus obtained, we compare the result with the one obtained with the help of publicly available global optimization software.
We implement the energy function formulation in IPOPT [23], which is an interior point method for local optimization. We provide IPOPT with first derivative information by using Lemma 1(a), and allow for automatic numerical Hessian computations. For the MIQP implementation, we use BONMIN [24], which is a mixed integer convex programming software that uses IPOPT as the local search subroutine. We use SCIP [25], a general purpose constraint optimization software, for global optimization.
Table I shows the cost comparison between the energy function formulation, the MIQP relaxation and SCIP. The solution provided by the energy function formulation is a feasible solution and hence it is an upper bound. We conclude, based on the experiments, that both the lower bound obtained by the MIQP and the upper bound from the energy function formulation are fairly close to the global optimum.
Pressure Bounds       

Energy function heuristics  651  576  491 
Global Solver  665  579  501 
MIQP  685  581  501 
Table II shows the optimal cost achieved by the energy function formulation, the MIQP relaxation and by SCIP. Similar to the case without compression, one observes that the lower and upper bounds found by the MIQP and by the energy function based formulation are close to the global optimum.
Pressure Bounds       

Energy function heuristics  694  627  584 
Global Solver  737  632  598 
MIQP  752  640  598 
Viii Conclusions and Path Forward
We developed a strictly convex variational representation of dissipative network flow equations that proves the uniqueness of the solution and provides an efficient algorithm to find it numerically. For the maximum throughput problem, we presented a twice differentiable optimization formulation using energy functions as a method to obtain a solution, and a mixed integer convex relaxation to test the optimality (or degree of suboptimality) of the solution. We test both methods on a model of natural gas networks and observe that they were successful in producing solutions close to global optimality.
For future work, it would be useful to extend these techniques to dynamic optimization problems on dissipative flow networks (such as transient gas networks). The optimization formulations could also be used as subroutines in multilevel programs associated with network planning and expansion. Coupled with the monotonicity property, it also provides scope for extension to network optimization problems accounting for uncertainty such as the robust maximum throughput problem [3] and the problem of finding the most probable fluctuations that cause violation of bounds known as the instanton problem (see e.g. a similar setting discussed in the context of power systems [26]).
Ix Acknowledgements
The authors thank S. Backhaus for multiple discussions and advice, and A. Zlotnik for providing examplary model of the gas network he designed. The work at LANL was carried out under the auspices of the National Nuclear Security Administration of the U.S. Department of Energy at Los Alamos National Laboratory under Contract No. DEAC5206NA25396 and it was partially supported by DTRA Basic Research Project . The authors also acknowledge partial support of the Advanced Grid Modeling Program in the US Department of Energy Office of Electricity.
Proof of Monotonicity Property (Corollary 1): By Lemma 1(a), , which means . By Lemma 1(b), we have the . From (23) and (25) we observe that the Hessian is a diagonally dominant invertible matrix with positive diagonal entries and nonpositive offdiagonal entries. From a standard result in linear algebra we have that all entrees in are nonnegative.
References
 [1] J. J. Maugis, “Etude de rÃ©seaux de transport et de distribution de fluide [in french],” pp. 243â–248, 1977.
 [2] F. Babonneau, Y. Nesterov, and J.P. Vial, “Design and operations of gas transmission networks,” Operations Research, 2012. [Online]. Available: http://or.journal.informs.org/content/early/2012/02/10/opre.1110.1001.abstract
 [3] M. Vuffray, S. Misra, and M. Chertkov, “Monotonicity of dissipative flow networks renders robust maximum profit problem tractable: General analysis and application to natural gas flows,” in Submitted to 54th IEEE Conference on Decision and Control, Osaka, 2015.
 [4] B. Bollobas, Modern Graph Theory. Springer, 1998.
 [5] P. Aylett, “The energyintegral criterion of transient stability limits of power systems,” Proceedings of the IEEPart C: Monographs, vol. 105, no. 8, pp. 527–536, 1958.
 [6] A. Araposthatis, S. Sastry, and P. Varaiya, “Analysis of powerflow equation,” International Journal of Electrical Power & Energy Systems, vol. 3, no. 3, pp. 115–126, 1981.
 [7] A. Bergen and D. Hill, “A structure preserving model for power system stability analysis,” Power Apparatus and Systems, IEEE Transactions on, vol. 100, no. 1, pp. 25–35, 1981.
 [8] K. Dvijotham, S. Low, and M. Chertkov, “Convexity of EnergyLike Functions: Theoretical Results and Applications to Power System Operations,” arxiv:1501.04052, 2015.
 [9] A. Osiadacz, Simulation and analysis of gas networks. Gulf Publishing Company,Houston, TX, Jan 1987.
 [10] S. Wu, R. RíosMercado, E. Boyd, and L. Scott, “Model relaxations for the fuel cost minimization of steadystate gas pipeline networks,” Mathematical and Computer Modelling, vol. 31, no. 23, pp. 197 – 220, 2000. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0895717799002320
 [11] S. Misra, M. Fisher, S. Backhaus, R. Bent, M. Chertkov, and F. Pan, “Optimal compression in natural gas networks: A geometric programming approach,” Control of Network Systems, IEEE Transactions on, vol. 2, no. 1, pp. 47–56, March 2015.
 [12] G. Como, K. Savla, D. Acemoglu, M. A. Dahleh, and E. Frazzoli, “On robustness analysis of largescale transportation networks,” in Proceedings of the International Symposium on Mathematical Theory of Networks and Systems, 2010, pp. 2399–2406.
 [13] P. Varaiya, “Max pressure control of a network of signalized intersections,” Transportation Research Part C: Emerging Technologies, vol. 36, pp. 177–195, 2013.
 [14] ——, “The maxpressure controller for arbitrary networks of signalized intersections,” in Advances in Dynamic Network Modeling in Complex Transportation Systems. Springer, 2013, pp. 27–66.
 [15] S. Boyd and L. Vandenberghe, Convex Optimization. Cambridge University Press, 2004.
 [16] R. T. Rockafellar, Convex Analysis. Princeton University Press, 1970.
 [17] R. Merris, “Laplacian matrices of graphs: a survey,” Linear Algebra Applications, vol. 197, 1994.
 [18] S. Chaiken and D. J. Kleitman, “Matrix tree theorems,” Journa of Combinatorial Theory, vol. 24, no. 377381, 1978.
 [19] G. McCormick, “Computability of global solutions to factorable nonconvex programs: Part I  convex underestimating problems,” Mathematical Programming, vol. 10, no. 146175, 1976.
 [20] P. Wong and R. Larson, “Optimization of naturalgas pipeline systems via dynamic programming,” Automatic Control, IEEE Transactions on, vol. 13, no. 5, pp. 475–481, 1968.
 [21] C. BorrazSánchez, “Optimization methods for pipeline transportation of natural gas,” Ph.D. dissertation, Department of Informatics, University of Bergen, Norway, October 2010.
 [22] S. Misra, M. W. Fisher, R. Bent, S. Backhaus, M. Chertkov, and F. Pan, “Optimal compression in natural gas networks: a geometric programming approach,” IEEE Control of Network Systems, 2014.
 [23] A. W. L. T. Biegler, “On the implementation of a primaldual interior point filter line search algorithm for largescale nonlinear programming,” Mathematical Programming, vol. 106, no. 1, pp. 25–27, 2006.
 [24] P. Bonami, A. Wachter, L. T. Biegler, A. R. Conn, G. Cornuejols, I. E. Grossman, C. D. Laird, J. Lee, A. Lodi, F. Margot, and N. Sawaya, “An algorithmic framework for convex mixed integer nonlinear programs.” IBM Research Report, RC23771, 2005.
 [25] T. Achterberg, T. Berthold, T. Koch, and K. Wolter, “Constraint integer programming: a new approach to integrate CP and MIP,” Integration of AI and OR techniques in constraint programming for combinatorial optimization problems, CPAIOP, 2008.
 [26] M. Chertkov, F. Pan, and M. Stepanov, “Predicting failures in power grids: The case of static overloads,” Smart Grid, IEEE Transactions on, vol. 2, no. 1, pp. 162–172, March 2011.