# Stationary Gas Networks with Compressor Control and Random Loads: Optimization with Probabilistic Constraints

Authors: Martin Gugat; Michael Schuster
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

We introduce a stationary model for gas flow based on simplified isothermal Euler equations in a non-cycled pipeline network. Especially the problem of the feasibility of a random load vector is analyzed. Feasibility in this context means the existence of a flow vector meeting these loads, which satisfies the physical conservation laws with box constraints for the pressure. An important aspect of the model is the support of compressor stations, which counteract the pressure loss caused by friction in the pipes. The network is assumed to have only one influx node; all other nodes are efflux nodes. With these assumptions the set of feasible loads can be characterized analytically. In addition we show the existence of optimal solutions for some optimization problems with probabilistic constraints. A numerical example based on real data completes this paper.


---

## Body

## 1. Introduction

In this paper, we present a way to model stationary gas networks with random loads. Natural gas as energy source has been popular for decades. Because of the nuclear power phase-out and aim to leave energy gained by coal behind, natural gas as energy source is more current than ever. So in this paper we study the problem of gas transport in a pipeline network mathematically. The aim of this paper is to get results in optimization problems with probabilistic constraints using the example of gas networks.There are several studies about the mathematical problem of gas transport. A gas transport network in general can be modeled as a system of hyperbolic balance laws, like e.g., the isothermal Euler equations. In [1] one can find a great overview about existing models and their application areas. A network system is modeled as a graph with nodes end edges. It can be extended by different elements, like compressor stations, valves and resistors. In [2] the authors describes the functionality of these elements and how they can be included to the model.We assume an active stationary state for our model that means we use compressor stations as controllable elements in a state of time-independent flows and pressures. Our model is based on the Weymouth-Equation (see [3]), a simplification of the isothermal Euler equation (see e.g., [4, 5]). Further [3] also studies the isothermal Euler equations for real gas. We assume our gas to be ideal. The optimal control problem in pipeline networks has been studied in, e.g., [4, 6]. In our work we assume the loads to be random, so this leads to optimization problems with probabilistic constraints (see [7]).In the next section, we introduce our gas network model, which is based on the model in [8]. This is a passive one which we extend by compressor stations. So our model contains inner control variables which makes our model an active one. These controls will be one goal of the optimization in Section 4.In Section3 we characterize the set of feasible loads. Instead of searching a flow and a pressure vector (based on a given load vector), which fulfills the physical balance laws, we look for these vectors fulfilling a system of inequalities, which depend on bounds for the pressure. In addition, we assume the load vector to be random. So we adapt the system to the spheric-radial decomposition, which is our central tool to handle the uncertainty in the loads.In Section4, we consider some optimization problems with and without probabilistic constraints. The existence of optimal solutions is based on the analysis we did in Section 3. This should built a base for further works about optimal control problems with probabilistic constraints using the example of gas networks.In Section5 we present three numerical examples. The first is a minimal graph without inner control, which shows the idea of the spheric-radial decomposition. The second one is a minimal graph with inner control and the third one is based on real data, for which we present real values for a network.In the last section we present a few ideas of extending this paper. One main aspect is the turnpike-theory.


## 2. Mathematical Modeling

We start with the introduction of the model that is also used in [8]. This model does not include compressor stations. We extend this model to include compressor stations and we discuss the network analysis. Then the main tool, the so called spheric-radial decomposition, is introduced to handle the stochastic uncertainty.


### 2.1. Model Description

We consider a connected, directed graphG=(V+,E) which represents a pipeline gas transport network. We assume that there are no cycles in the graph, so the network is a tree. We set V+≔n+1 and E≔m=n (n,m∈N). In our model, every node is either a influx node (gas enters the network) or an efflux node (gas leaves the network). An edge can either be a flux edge, so the pressure decreases along the edge, or a compressor edge, so the pressure increases along the edge. We define EF as the set of all flux edges and EC as the set of compressor edges, it follows E=EF∪EC with EF∩EC=∅.Letb+∈Rn+1 with 1n+1Tb+=0 denote the balanced load vector and assume bi≤0 for nodes with gas influx and bi≥0 for nodes with gas efflux (i∈{0,…,n}). The vector 1n+1 denotes the vector of all ones in the dimension n+1. Furthermore, a vector q∈Rm describes the flows through the edges resp. through the pipes. The pressures at the nodes are defined in a vector p+∈Rn+1. Let pressure bounds p+,min,p+,max∈Rn+1 be given. For the pressure we consider the constraints p+∈[p+,min,p+,max]. For further modeling we need the following definition:Definition 1.
Consider the graphG=(V+,E):(i)
h(e) denotes the head node of an edge e and f(e) denotes the feet of an edge e for all e∈E.(ii)
E0(v)≔{e∈E∣h(e)=v∨f(e)=v} denotes all edges which are connected to node v.(iii)
The matrixAij+∈Rn+1×m, Ai,j+≔σ(vi,ej) with(1)σv,e≔-1iffe=vande∈E0v1ifhe=vande∈E0v0ife∉E0vis called the incidence matrix of the graphG.(iv)
Forv∈V+ is Π(v)≔{v}∪⋃{k∈V+∣∃e∈E:f(e)=k∧h(e)=v}Π(k) the (unique) directed path from the root to v.(v)
Forv,w∈V+ is Πw(v)≔Π(v)∖Π(w) the (unique) directed path from w to v.Note that definition (iv) only makes sense since we assumed the graph to be tree-structured, so there are no cycles and the union is finite. The model is based upon the conservation equations of mass and momentum. The mass equation for the nodes is formulated for every node by(2)∑e∈E0vqe=bv+∀v∈V+.This is equivalent toKirchhoff’s first law (see [9]) and by using the definition of the incidence matrix, the equation for mass conservation for the whole graph is(3)A+q=b+resp.Aq=bwithA+=A0A∈R1×m∈Rn×mandb+=b0b∈R∈Rn.As balance laws we use the isothermal Euler equations (see [3, 4, 10]) for a horizontal pipe, so the flow through every flux edge is modeled by(4)ρte+qxe=0qte+qe2ρe+a2ρex=-λe2Deqeqeρe∀e∈EFHere, ρ is the gas density, a is the speed of sound, λ is a (constant) friction coefficient and D is the diameter of the pipe.We assume first, that the network is in a steady state, so the time derivatives are equal to zero. And second, we assume the gas flow to be slow, so the coefficient(qe/a)2 is negligible. These assumptions simplify the equations enormously. With the ideal gas equation p=ρRST (see [9]) we get a direct dependency between pressure and density. Solving this simplified equation leads us to the so called Weymouth-Equation (see [3]):(5)pfe2-phe2=ϕeqeqewithϕe=RST2λeDeLea2∀e∈EF.Here, RS is the specific gas constant, T is the temperature of the gas and Le is the length of the pipe e. This equation shows the pressure drop along the pipes (see Figure 1). Some articles consider the isothermal Euler equations in gas transport without these simplifications, e.g., [5, 10].Figure 1
Pressure drop along a pipe forD=0.762m, T=289.15K, RS=500.27J/kgK, λ=0.005, a=429.17m/s, q=795kg/m2s and p(0)=60⋅105Pa.Compressor stations counteract the pressure drop along the pipes. These stations are modeled as pipes without friction. We use the following model equation (see [1], section 7.4.7; [11], section 2.3):(6)cγpoutpinγ-1=u~.Here, pout and pin are the pressures at the end resp. at the beginning of one pipe, c and γ are constants and u~, which is our control, is the specific change in adiabatic enthalpy. More information about the parameters can be found in [11]. In addition, [2] contains more information about the functionality of compressor stations and a deduction of the equation. A short transformation of the equation leads to our representation of the pressure rise along the compressor edges:(7)phepfe2=uewithue=γcu~+12/γ∀e∈EC.It holds ue≥1 where ue=1 means that the compressor station is switched off (see Figure 2).Figure 2
Compressor property forc=144652.6J/kg, γ=0.225, and pin=40⋅105Pa.Now, we define the set of feasible loads. We say that a vectorb+∈Rn+1 is feasible, if (3), (5), and (7) are fulfilled. So the feasible set (here called M^) is defined as follows:(8)M^=b+∈Rn+1∣1n+1Tb+=0and∃p+,q∈Rn+1×Rm:p+∈p+,min,p+,maxand3,5,7are  fulfilled.


### 2.2. Stochastic Uncertainty in the Loads

In reality, the future demand for gas can never be known exactly a priori. It depends on various physical effects like e.g., the temperature of the environment. Also personal reasons can be responsible for a higher or lower use of gas. Thus the load vector is considerd to be a random vector. Many papers study how the distribution functions of the load vector can be identified, e.g., [12]. This article also describes, why a multivariate Gaussian distribution is a good choice for random gas demand. So we use some random vector ξ(ω) with ξ~N(μ,Σ) for a mean vector μ and a positive definite covariance matrix Σ on a suitable probability space (Ω,A,P) and our aim is to calculate the following probability:(9)Pω∈Ω∣ξω∈Mfor a suitable set M.In our case, we assume that the tree-structured graph has only one influx node, which gets the number zero, all remaining nodes are efflux nodes. Then the graph is numerated with breadth-first search or depth-first search and every edgee∈E gets the number max{h(e),f(e)}. Our set M is defined as follows:(10)M≔b∈R≥0n∣-1nTb,b∈M^.There are many studies for special setsM, e.g., for polyhedral sets (see [13]), but our set is very general. So we need a very general method for calculating the probability (9). Our method of choice is the so calledspheric-radial decomposition. It became popular for dealing with probabilistic constraints (e.g., see [8, 14–19]).Theorem 2 (spheric-radial decomposition).
Letξ~N(0,R) be the n-dimensional standard Gaussian distribution with zero mean and positive definite correlation matrix R. Then, for any Borel measurable subset M⊆Rn it holds that(11)Pξ∈M=∫Sn-1μχr≥0∣rLv∈Mdμηv,where Sn-1 is the (n-1)-dimensional sphere in Rn, μη is the uniform distribution on Sn-1, μχ denotes the χ-distribution with n degrees of freedom and L is such that R=LLT (e.g., Cholesky decomposition).The difference with our problem is, that our load vector is normal distributed with mean valueμ and positive definite covariance Σ. We want to use the result of Theorem 2 to general Gaussian distributions (see [8]). Therefore we set(12)ξ∗≔D-1ξ-μ~N0,RwithD≔diagΣiii=1,…,nandR≔D-1ΣD-1,and observe that P(ξ∈M)=P(ξ∗∈D-1(M-μ)). So our aim is to compute the following probability:(13)Pξ∗∈D-1M-μ.So the calculation of the probability (9) for a very general set M is equivalent to compute the integral in Theorem 2, which is not easy to compute analytically, too. But the integral can be computed numerically in a easy and nearly exact way. Therefore we formulate the following algorithm (see [8], Algorithm 4):Algorithm 3.
Letξ~N(μ,Σ) and L such that LLT=Σ.(1)
SampleN points {v1,…,vN} uniformly distributed on the sphere Sn-1.(2)
Compute the one-dimensional setsMi≔{r≥0∣rLvi+μ∈M} for i=1,…,N.(3)
SetP(ξ∈M)≈1/N∑i=1Nμχ(Mi).The first step is quite easy. There are many possibilities to sample a set of uniformly distributed points on the sphereSn-1, e.g.,Monte-Carlo sampling or therandn-generator inMATLAB®. Normalizing these points leads us to the sampling we are looking for. Alternatively, quasi Monte-Carlo methods can be used. In [8] the authors use a quasi Monte-Carlo method for their computations and they show its advantages in Figure 6.The main challenge is to handle the one-dimensional sets of the second step. If the setM is convex, these sets are closed intervals but M need not to be convex. The aim is to get a representation of Mi (i=1,…,N) as a union of disjoint intervals. With this, the third step can be done efficiently. The next section deals with this representation in detail.The third step is the approximation of the integral. The valuesμχ(Mi) can be calculated easily with the distribution function Fχ of the distribution (see [20] for details of that distribution). With the representation Mi=⋃j=1t[a_j,a¯j] the χ-distribution can be calculated the following:(14)μχMi=∑j=1tFχa¯j-Fχa_j.This leads to the wanted probability. The χ-distribution function is given by the following (see [20], p. 132):(15)Fχx=12N/2ΓN/2∫0xvN/2-1e-v/2dvwith Γ(1/2)=π, Γ(1)=1 and Γ(r+1)=rΓ(r) for all r∈R+. Note that there are very precise numerical approximations for the distribution function, e.g., thechi2cdf-function inMATLAB® (see Figure 3).Figure 3
Distribution function of the Chi-square distribution for 3 (blue), 5 (red), 8 (green), and 12 (purple) degrees of freedom.


## 2.1. Model Description

We consider a connected, directed graphG=(V+,E) which represents a pipeline gas transport network. We assume that there are no cycles in the graph, so the network is a tree. We set V+≔n+1 and E≔m=n (n,m∈N). In our model, every node is either a influx node (gas enters the network) or an efflux node (gas leaves the network). An edge can either be a flux edge, so the pressure decreases along the edge, or a compressor edge, so the pressure increases along the edge. We define EF as the set of all flux edges and EC as the set of compressor edges, it follows E=EF∪EC with EF∩EC=∅.Letb+∈Rn+1 with 1n+1Tb+=0 denote the balanced load vector and assume bi≤0 for nodes with gas influx and bi≥0 for nodes with gas efflux (i∈{0,…,n}). The vector 1n+1 denotes the vector of all ones in the dimension n+1. Furthermore, a vector q∈Rm describes the flows through the edges resp. through the pipes. The pressures at the nodes are defined in a vector p+∈Rn+1. Let pressure bounds p+,min,p+,max∈Rn+1 be given. For the pressure we consider the constraints p+∈[p+,min,p+,max]. For further modeling we need the following definition:Definition 1.
Consider the graphG=(V+,E):(i)
h(e) denotes the head node of an edge e and f(e) denotes the feet of an edge e for all e∈E.(ii)
E0(v)≔{e∈E∣h(e)=v∨f(e)=v} denotes all edges which are connected to node v.(iii)
The matrixAij+∈Rn+1×m, Ai,j+≔σ(vi,ej) with(1)σv,e≔-1iffe=vande∈E0v1ifhe=vande∈E0v0ife∉E0vis called the incidence matrix of the graphG.(iv)
Forv∈V+ is Π(v)≔{v}∪⋃{k∈V+∣∃e∈E:f(e)=k∧h(e)=v}Π(k) the (unique) directed path from the root to v.(v)
Forv,w∈V+ is Πw(v)≔Π(v)∖Π(w) the (unique) directed path from w to v.Note that definition (iv) only makes sense since we assumed the graph to be tree-structured, so there are no cycles and the union is finite. The model is based upon the conservation equations of mass and momentum. The mass equation for the nodes is formulated for every node by(2)∑e∈E0vqe=bv+∀v∈V+.This is equivalent toKirchhoff’s first law (see [9]) and by using the definition of the incidence matrix, the equation for mass conservation for the whole graph is(3)A+q=b+resp.Aq=bwithA+=A0A∈R1×m∈Rn×mandb+=b0b∈R∈Rn.As balance laws we use the isothermal Euler equations (see [3, 4, 10]) for a horizontal pipe, so the flow through every flux edge is modeled by(4)ρte+qxe=0qte+qe2ρe+a2ρex=-λe2Deqeqeρe∀e∈EFHere, ρ is the gas density, a is the speed of sound, λ is a (constant) friction coefficient and D is the diameter of the pipe.We assume first, that the network is in a steady state, so the time derivatives are equal to zero. And second, we assume the gas flow to be slow, so the coefficient(qe/a)2 is negligible. These assumptions simplify the equations enormously. With the ideal gas equation p=ρRST (see [9]) we get a direct dependency between pressure and density. Solving this simplified equation leads us to the so called Weymouth-Equation (see [3]):(5)pfe2-phe2=ϕeqeqewithϕe=RST2λeDeLea2∀e∈EF.Here, RS is the specific gas constant, T is the temperature of the gas and Le is the length of the pipe e. This equation shows the pressure drop along the pipes (see Figure 1). Some articles consider the isothermal Euler equations in gas transport without these simplifications, e.g., [5, 10].Figure 1
Pressure drop along a pipe forD=0.762m, T=289.15K, RS=500.27J/kgK, λ=0.005, a=429.17m/s, q=795kg/m2s and p(0)=60⋅105Pa.Compressor stations counteract the pressure drop along the pipes. These stations are modeled as pipes without friction. We use the following model equation (see [1], section 7.4.7; [11], section 2.3):(6)cγpoutpinγ-1=u~.Here, pout and pin are the pressures at the end resp. at the beginning of one pipe, c and γ are constants and u~, which is our control, is the specific change in adiabatic enthalpy. More information about the parameters can be found in [11]. In addition, [2] contains more information about the functionality of compressor stations and a deduction of the equation. A short transformation of the equation leads to our representation of the pressure rise along the compressor edges:(7)phepfe2=uewithue=γcu~+12/γ∀e∈EC.It holds ue≥1 where ue=1 means that the compressor station is switched off (see Figure 2).Figure 2
Compressor property forc=144652.6J/kg, γ=0.225, and pin=40⋅105Pa.Now, we define the set of feasible loads. We say that a vectorb+∈Rn+1 is feasible, if (3), (5), and (7) are fulfilled. So the feasible set (here called M^) is defined as follows:(8)M^=b+∈Rn+1∣1n+1Tb+=0and∃p+,q∈Rn+1×Rm:p+∈p+,min,p+,maxand3,5,7are  fulfilled.


## 2.2. Stochastic Uncertainty in the Loads

In reality, the future demand for gas can never be known exactly a priori. It depends on various physical effects like e.g., the temperature of the environment. Also personal reasons can be responsible for a higher or lower use of gas. Thus the load vector is considerd to be a random vector. Many papers study how the distribution functions of the load vector can be identified, e.g., [12]. This article also describes, why a multivariate Gaussian distribution is a good choice for random gas demand. So we use some random vector ξ(ω) with ξ~N(μ,Σ) for a mean vector μ and a positive definite covariance matrix Σ on a suitable probability space (Ω,A,P) and our aim is to calculate the following probability:(9)Pω∈Ω∣ξω∈Mfor a suitable set M.In our case, we assume that the tree-structured graph has only one influx node, which gets the number zero, all remaining nodes are efflux nodes. Then the graph is numerated with breadth-first search or depth-first search and every edgee∈E gets the number max{h(e),f(e)}. Our set M is defined as follows:(10)M≔b∈R≥0n∣-1nTb,b∈M^.There are many studies for special setsM, e.g., for polyhedral sets (see [13]), but our set is very general. So we need a very general method for calculating the probability (9). Our method of choice is the so calledspheric-radial decomposition. It became popular for dealing with probabilistic constraints (e.g., see [8, 14–19]).Theorem 2 (spheric-radial decomposition).
Letξ~N(0,R) be the n-dimensional standard Gaussian distribution with zero mean and positive definite correlation matrix R. Then, for any Borel measurable subset M⊆Rn it holds that(11)Pξ∈M=∫Sn-1μχr≥0∣rLv∈Mdμηv,where Sn-1 is the (n-1)-dimensional sphere in Rn, μη is the uniform distribution on Sn-1, μχ denotes the χ-distribution with n degrees of freedom and L is such that R=LLT (e.g., Cholesky decomposition).The difference with our problem is, that our load vector is normal distributed with mean valueμ and positive definite covariance Σ. We want to use the result of Theorem 2 to general Gaussian distributions (see [8]). Therefore we set(12)ξ∗≔D-1ξ-μ~N0,RwithD≔diagΣiii=1,…,nandR≔D-1ΣD-1,and observe that P(ξ∈M)=P(ξ∗∈D-1(M-μ)). So our aim is to compute the following probability:(13)Pξ∗∈D-1M-μ.So the calculation of the probability (9) for a very general set M is equivalent to compute the integral in Theorem 2, which is not easy to compute analytically, too. But the integral can be computed numerically in a easy and nearly exact way. Therefore we formulate the following algorithm (see [8], Algorithm 4):Algorithm 3.
Letξ~N(μ,Σ) and L such that LLT=Σ.(1)
SampleN points {v1,…,vN} uniformly distributed on the sphere Sn-1.(2)
Compute the one-dimensional setsMi≔{r≥0∣rLvi+μ∈M} for i=1,…,N.(3)
SetP(ξ∈M)≈1/N∑i=1Nμχ(Mi).The first step is quite easy. There are many possibilities to sample a set of uniformly distributed points on the sphereSn-1, e.g.,Monte-Carlo sampling or therandn-generator inMATLAB®. Normalizing these points leads us to the sampling we are looking for. Alternatively, quasi Monte-Carlo methods can be used. In [8] the authors use a quasi Monte-Carlo method for their computations and they show its advantages in Figure 6.The main challenge is to handle the one-dimensional sets of the second step. If the setM is convex, these sets are closed intervals but M need not to be convex. The aim is to get a representation of Mi (i=1,…,N) as a union of disjoint intervals. With this, the third step can be done efficiently. The next section deals with this representation in detail.The third step is the approximation of the integral. The valuesμχ(Mi) can be calculated easily with the distribution function Fχ of the distribution (see [20] for details of that distribution). With the representation Mi=⋃j=1t[a_j,a¯j] the χ-distribution can be calculated the following:(14)μχMi=∑j=1tFχa¯j-Fχa_j.This leads to the wanted probability. The χ-distribution function is given by the following (see [20], p. 132):(15)Fχx=12N/2ΓN/2∫0xvN/2-1e-v/2dvwith Γ(1/2)=π, Γ(1)=1 and Γ(r+1)=rΓ(r) for all r∈R+. Note that there are very precise numerical approximations for the distribution function, e.g., thechi2cdf-function inMATLAB® (see Figure 3).Figure 3
Distribution function of the Chi-square distribution for 3 (blue), 5 (red), 8 (green), and 12 (purple) degrees of freedom.


## 3. Explicit Characterization of the Gas Flow

In this section we want to characterize the feasible set to use Algorithm3. The aim is to get a representation of the one-dimensional sets of step (2) as a union of disjoint intervals. Therefore we need another representation of the feasible set as a system of inequalities. This is an idea from [8], Theorem 1. In the proof, the authors use a representation of the pressure loss for the full network graph:(16)A+Tp+2=-Φqq,which is equivalent to (5) if and only if EC=∅. In this paper we also consider the case EC≠∅. The matrix Φ is diagonal with the values ϕe as defined in (5) at its diagonal. This representation is not possible yet, because we extended our model by compressor edges, so this equation does not hold for all edges. So we have to find a way to use a similar equation like (16) to state a suitable theorem. The idea that allows to include compressor stations is to remove compressor edges from the graph and replace them by suitable algebraic conditions. First we consider the equation of mass:(17)∑e∈E0vqe=bv+⇔∑e∈E0v∩EFqe=bv+-∑e∈E0v∩ECqe.We want to treat the flows through the compressor edges as loads, so we have to adapt the load vector. The new load vector b~+ is defined as follows:(18)b~v+≔bv++qe∀e∈E0vwithe∈ECandfe=vbv+-qe∀e∈E0vwithe∈ECandhe=vbv+else.The resulting graph G~=(V+,EF) is no longer connected (scheme is shown in Figure 4).Figure 4
Network graph with two compressor edges (left) and new not connected network graph with three connected subgraphs after removing the compressor edges (right).We setEF≔m1 and EC≔m2 with m1+m2=m, so G~ consists of m2+1 subgraphs, which are all connected. Note that q is a vector in Rm1 now. We determine that Gi is the i-th self-connected subgraph of G (depends on numbering of the graph). Analogously pi,j+ is the j-th component of the i-th subgraph. The terms qi,j and Ai,j,k are defined similarly. The feasible set for G~ is then(19)M~=b~+∈Rn+1∣1n+1Tb^+=0∧∃p+,q∈R≥0n+1×Rmwithp+∈p+,min,p+,max:Aiqi=b~i∀i=1,…,m2+1Ai+Tpi+2=-Φiqi∘qi∀i=1,…,m2+1phepfe2=ue∀e∈ECThe last equation inM~ is the same as in (8), that models the pressure rise caused by the compressor. With this we can present the equation for pressure loss in the form of (16) for every subgraph. Before we state a theorem for characterizing the feasible set, we show the equivalence of M^ and M~.Lemma 4.
Letb~+ be as defined in (18). The feasible sets defined in (8) and (19) are equivalent in the following sense:(20)iFor  allb+∈M^there  isb~+∈M~iiFor  allb~+∈M~there  isb+∈M^Proof.
(i) Consider a vectorb+∈M^. It holds 1n+1Tb+=0 and because of (18) it follows directly 1n+1Tb~+=0. After removing the compressor edges, from A+q=b+ it follows with (17) that Ai+qi=b~i+, which is equivalent to Aiqi=b~i. The pressure loss is fulfilled for all flux edges in M^ and because the flow vector q doesn’t change in the relevant coefficients while removing the compressor edges, the pressure loss is fulfilled in M~ too. The compressor property is independent of b+ and q, so it also holds. From this it follows b~+∈M~.
(ii) The other way around is quite similar. We consider a vectorb~+∈M~. From 1n+1Tb~+=0 it follows with (18) that 1n+1Tb+=0. The equation for mass conservation is fulfilled in every subgraph. So with (17) it follows A+q=b+, which implies the mass conservation in M. The pressure loss holds in all self-connected parts of the graph too and because the relevant coefficients of q doesn’t change while adding the compressor edges, the pressure loss is fulfilled in M too. Last, as said before, the compressor property is independent of b~+ and q, it is also fulfilled in M. This completes the proof.Because our network graph is a tree, the mass conservation is fulfilled with a full rank incidence matrixA of size n×n, so q can be computed explicit by the following:(21)q=A-1b,which makes us able to use Lemma 4.With this result, we can write the equation for momentum conservation depending on the incidence matrices for the parts of the tree. This makes us able to formulate the theorem for characterizing the feasible setM~. Therefore we will use the following function:(22)g:Rn→Rx↦AT-1ΦA-1xA-1x,which describes the pressure loss between the nodes. We also introduce a new notation:(i)
u(i,j),k is the k-th control on the path from vi,0 to vj,0(ii)
g(i,j),k is the function g for the k-th subgraph on the path from vi,0 to vj,0, whereas g(i,j),1 is the function g of the graph Gi. According to this, g(i,j),k,l is the l-th component of the function g(i,j),k(iii)
b~(i,j),k is the load vector of the k-th subgraph in the path from vi,0 to vj,0, whereas b~(i,j),1 is the load vector of the graph Gi and b~(i,j),k,l is the l-th component of b~(i,k),kAnalogously we defineA(i,j),k, p(i,j),k and q(i,j),k. Furthermore we define the following values:(i)
ki,j∗≔max{k∈{1,…,m2+1}∣vk,0∈Π(vi,0)∧vk,0∈Π(vj,0)} as the largest index of all subgraph, the paths to vi,0 and vj,0 pass through(ii)
ni,j∗:k∈1,…,m2+1∣vk,0∈Πvi,0vj,0 as the number of subgraph which are on the path from vi,0 from vj,0 containing Gi and Gj(iii)
mi,j∗≔e∈EC∣he∈Πvi,0vj,0∧fe∈Πvi,0vj,0 as the number of controls which are between vi,0 and vj,0With this notation and these values we define a sum fori,j∈{1,m2+1}. If vi,0∈Π(vj,0) we set ∑i,i≔0 and else we set(23)Σk∗,i≔∑k=1n∗-21∏l=1m∗-kuk∗,i,lgk∗,i,n∗-k,feuk∗,i,n∗-kb~k∗,i,n∗-k+gk∗,i,1,feuk∗,i,1b~k∗,i,1,where we write only n∗ and m∗ instead of nk∗,i∗ and mk∗,i∗ for better reading. Last, we define a product Πi,i≔1 and(24)Πk∗,i≔∏k=1m∗uk∗,i,k.From a first point of view, the notation seems to be confusing, but using this notation makes sense since we want to guarantee the feasibility of a load vector by comparing the pressure bounds at every node with each other. We explain this in a short example. Figure 5 shows a graph after removing the compressor edges.Figure 5
Example graph for illustration of the used notation.Figure 6
Graph of Example1.If we want to compare the pressure bounds of node 11 and node 12, we need to know how the pressures change on the path from 2 to 11, resp., from 2 to 12. Node 11 is part of the subgraphG5 and node 12 is part of the subgraph G6, so it follows k5,6∗=2 because node 2 is part of subgraph 2. The notation can be understood as a numbering along the path. It is u(2,11),1=u2, the first control on the path from node 2 to node 11 and u(2,11),2=u4, the second control on this path. Further it follows n2,11∗=3 (G2,G3,G5) and m2,11∗=2 (u2,u4). The product defined before is the product of all controls along this path and the sum is the pressure loss between node 2 and node 7. Using this notation makes us able to formulate the theorem for characterizing the feasible set in a readable way.Theorem 5.
For given pressure boundsp+,min,p+,max∈Rn+1 and controls ui∈R (i=1,…,m2) the following equivalence holds.
A vectorb~+ with 1n+1Tb~+=0 is feasible, i.e., b~+∈M~, if and only if the following inequalities hold. For all i=1,…,m2+1 holds (feasibility inside the subgraphs):(25)pi,0min2≤mink=1,…,ni⁡pi,kmax2+gi,kb~i,(26)pi,0max2≥maxk=1,…,nipi,kmin2+gi,kb~i,(27)maxk=1,…,nipi,kmin2+gi,kb~i≤mink=1,…,nipi,kmax2+gi,kb~i.For all i,j=1,…,m2+1 with i<j holds (feasibility between the subgraphs):(28)1Πk∗,ipi,0min2+Σk∗,i≤1Πk∗,jpj,0max2+Σk∗,j,(29)1Πk∗,ipi,0max2+Σk∗,i≥1Πk∗,jpj,0min2+Σk∗,j,(30)1Πk∗,ipi,0min2+Σk∗,i≤1Πk∗,jmink=1,…,njpj,kmax2+gj,kb~j+Σk∗,j,(31)1Πk∗,ipi,0max2+Σk∗,j≥1Πk∗,jmaxk=1,…,njpj,kmin2+gj,kb~j+Σk∗,j,(32)1Πk∗,imaxk=1,…,nipi,kmin2+gi,kb~i+Σk∗,i≤1Πk∗,jpj,0max2+Σk∗,j,(33)1Πk∗,imink=1,…,nipi,kmax2+gi,kb~i+Σk∗,i≥1Πk∗,jpj,0min2+Σk∗,j,(34)1Πk∗,imaxk=1,…,nipi,kmin2+gi,kb~i+Σk∗,i≤1Πk∗,jmink=1,…,njpj,kmax2+gj,kb~j+Σk∗,j,(35)1Πk∗,imink=1,…,nipi,kmax2+gi,kb~i+Σk∗,i≥1Πk∗,jmaxk=1,…,njpj,kmin2+gj,kb~j+Σk∗,j.Let us have a look at the inequalities. Inequalities (25)–(27) are well-known from [8] and they describe the feasibility inside a connected tree. The other inequalities (28)–(35) guarantee the feasibility between the subgraphs. Therefore we compare the pressures of the entries of two parts ((28), (29)), of the entry of one part and the exit nodes of the other parts ((30)–(33)) and we compare the pressures of all exits ((34), (35)). Because we cannot compare the pressures of Gi and Gj directly (there may be no direct connection), we have to compare them using the subgraph Gk∗ which is the last part the paths to vi,0 and vj,0 pass. Therefore we defined the sums Σk∗,i and Σk∗,j, which trace back the pressures to Gk∗.Remark 6.
For an implementation it is wise to define a setU(i,j) containing all controls on the path from vi,0 to vj,0 and use an index running over this set to define the sums and the products.Now we want to proof the theorem.Proof of Theorem5.
Because the columns ofA+ sum up to zero, it holds(36)-1niTAi=Ai,0+i=1,…,m2+1.“⇒” We consider a feasible vector b+∈M~. From the equation of mass conservation of every subgraph it follows(37)qi=Ai-1b~i.We take the equation of momentum conservation for the i-th subgraph (i∈{1,…,m2+1}) and separate Ai+ to Ai,0+ and Ai:(38)Ai,0+Tpi,02+AiTpi2=-Φqiqi.Next we insert (36) into (38) and get(39)-AiT1nipi,02+AiTpi2=-Φqiqi.We multiply this equation from left with (AiT)-1 and insert (37), so it follows(40)1nipi,02=pi2+AiT-1ΦAi-1b~iAi-1b~iand using the definition of the function g we have(41)1nipi,02=pi2+gib~i.Note here, that pi is the vector of pressures of the subgraph Gi and gi(b~i) is the vector-valued function, which describes the pressure loss in Gi. Equation (41) holds for all subgraphs Gi with i∈{1,…,m2+1}. From this, the inequalities (25)–(27) follow directly with p+∈[p+,min,p+,max]. Next we use (41) component-by-component. For fix i,j∈{1,…,m2+1} with i<j it follows(42)pk∗,02=pk∗,k2+gk∗,kb~k∗k=1,…,nk∗.Remember, that ki,j∗ is the largest index of all subgraphs, the paths to vi,0 and vj,0 pass through. For better reading, we only write k∗ instead of ki,j∗. We use the equation of the compressor property for the first compressor on the path from vk∗,0 to vi,0,(43)pheuk∗,i,12pfeuk∗,i,12=uk∗,i,1,and insert it into (42) for k=f(eu(k∗,i),1):(44)pk∗,i,1,02=1uk∗,i,1pk∗,i,2,02+gk∗,i,1,feuk∗,i,1b~k∗,i,1.Note that pk∗,02=p(k∗,i),1,02, gk∗,k(b~k∗)=g(k∗,i),1,k(b~(k∗,i),1) and ph(eu(k∗,i),1)2=p(k∗,i),2,02. We use the compressor equation for the second compressor on the path from vk∗,0 to vi,0 and put it, together with (42) into (44):(45)pk∗,i,1,02=1uk∗,i,11uk∗,i,2pk∗,i,3,02+gk∗,i,2,feuk∗,i,2b~k∗,i,2+gk∗,i,1,feuk∗,i,1b~k∗,i,1.We repeat this procedure for all controls on the path from vk∗,0 to vi,0. Although we write n∗ and m∗ instead of nk∗,i∗ and mk∗,i∗. Thus we get(46)pk∗,i,1,02=1uk∗,i,11uk∗,i,2⋯1uk∗,i,m∗pk∗,i,n∗,02+gk∗,i,n∗-1,feuk∗,i,n∗-1b~k∗,i,n∗-1⋯+gk∗,i,2,feuk∗,i,2b~k∗,i,2+gk∗,i,1,feuk∗,i,1b~k∗,i,1.We take again (42) for the node vk∗,0 and do the same procedure just like before for the path from vk∗,0 to vj,0. With n∗=nk∗,j∗ and m∗=mk∗,j∗ we get(47)pk∗,j,1,02=1uk∗,j,11uk∗,j,2⋯1uk∗,j,m∗pk∗,j,n∗,02+gk∗,j,n∗-1,feuk∗,j,n∗-1b~k∗,j,n∗-1⋯+gk∗,j,2,feuk∗,j,2b~k∗,j,2+gk∗,j,1,feuk∗,j,1b~k∗,j,1.Expanding (46) and (47) and equalizing them (this is possible because p(k∗,i),1,0=p(k∗,j),1,0) leads to(48)1Πk∗,ipi,02+Σk∗,i=1Πk∗,jpj,02+Σk∗,jwith Σk∗,i resp. Σk∗,j defined before. (48) directly implies (28) and (29). Using (48) together with (42) for vi,0 first time, for vj,0 second time and for vi,0 and vj,0 a third time, the inequalities (30)–(35) follow directly. So this part of the proof is complete.
“⇐” For this part, we consider a vector b+∈Rn+1 with 1n+1Tb+=0, which fulfills the inequalities (25)–(35). We define the following sets (i=2,…,m2+1):(49)P1,1≔p1,0min2,p1,0max2,P2,1≔maxk=1,…,n1p1,kmin2+g1,kb~1,mink=1,…,n1p1,kmax2+g1,kb~1,P1,i≔1Π1,ipi,0min2+Σ1,i,1Πi,1pi,0max2+Σ1,i,P2,i≔1Π1,imaxk=1,…,nipi,kmin2+gi,kb~i+Σ1,i,1Πi,1mink=1,…,nipi,kmax2+gi,kb~i+Σ1,i.We set(50)P0≔⋂k=1m2+1P1,k∩⋂k=1m2+1P2,kand we will define a value p1,02∈P0 later. So to guarantee such a value p1,02 exists, we have to show first, that P0 is not empty. Because P0 is a finite intersection of convex intervals, we have to make sure that all intervals are not empty and every intersection between two intervals is not empty. In one dimension, this is sufficient that P0 is not empty. For i=1,…,m2+1, (51)P1,i≠∅accordingtorequirementP2,i≠∅becauseof26.Furthermore its P1,i∩P2,i≠∅ because of (25) and (26). The missing intersections are nonempty too, for i,j∈{1,…,m2+1} with i<j, as follows:(52)P1,i∩P1,j≠∅becauseof28and29.P1,i∩P2,j≠∅becauseof30and31.P2,i∩P1,j≠∅becauseof32and33.P2,i∩P2,j≠∅becauseof34and35.
So the setP0 is not empty. We define the following values:(53)p1,02∈P0,pi,02≔Π1,ip1,02-Σ1,ii=2,…,m2+1,pi,k2≔pi,02-gi,kb~ii=1,…,m2+1,k=1,…,ni,qi≔Ai-1b~ii=1,…,m2+1.Now we show that our choice of pressures and flows is feasible. It is p1,0∈[p1,0min,p1,0max] because of p1,02∈P1,1. Further because of p1,02∈P1,i for all i=2,…,m2+1 it holds(54)1Π1,ipi,02+Σ1,i≤1Π1,ipi,0max2+Σ1,iand(55)1Π1,ipi,02+Σ1,i≥1Π1,ipi,0min2+Σ1,i.Now we replace pi,02 in its definition with (pi,k2+gi,k(b~i)) and because of pi,02∈P2,i for all i=1,…,m2+1 it follows(56)1Π1,ipi,k2+gi,kb~i+Σ1,i≤1Π1,imink=1,…,nipi,kmax2+gi,kb~i+Σ1,iand(57)1Π1,ipi,k2+gi,kb~i+Σ1,i≥1Π1,imaxk=1,…,nipi,kmin2+gi,kb~i+Σ1,i.So for all i=1,…,m2+1 and for all k=1,…,ni it holds pi,k∈[pi,kmin,pi,kmax] which is equivalent to(58)p+∈p+,min,p+,max.The equation for mass conservation follows directly from the definition of qi (i=1,…,m2+1). For the pressure loss we use(59)pi2=1nipi,02-gib~i,And multiply this equation from left with AiT and use (36). With this it follows the equation for momentum conservation:(60)Ai+Tpi+2=-Φqiqi.The last equation missing is the compressor property. For a control j∈{1,…,m2} we fix a vertex vi,0 (i∈{1,…,m2+1}) with h(euj)=vi,0. We use pi,02=Π1,i(p1,02-Σ1,i) and solve the formula for p1,02 in the following way:(61)p1,02=1u1,i,1··1u1,i,m∗p1,i,n∗,02+g1,i,n∗-1,feu1,i,n∗-1b~1,i,n∗-1··+g1,i,1,feu1,i,1b~1,i,1.Consider w.l.o.g. that the subgraphs along the path from v1,0 to vi,0 are numbered chronological, s.t. if h(euj)∈Gi then f(euj)∈Gi-1. We use pi-1,02=Π1,i-1(p1,02-Σ1,i-1) to get the same representation for p1,02 depending on p(1,i-1),n∗,02, which is equivalent to p(1,i),n∗-1,0. We equalize these representations depending on pi,02 and pi-1,02 and we get(62)1u1,i,m∗p1,i,n∗,02+g1,i,n∗-1,feu1,i,n∗-1b~1,i,n∗-1=p1,i,n∗-1,02.We chose i∈{1,…,m2+1} s.t. uj is the last control on the path from v1,0 to vi,0, so it is u(1,i),m∗=uj. Further we have(63)p1,i,n∗,02=pi,02,p1,i,n∗-1,02=pi-1,02,g1,i,n∗-1,feu1,i,n∗-1b1,i,n∗-1=gi-1,feujbi-1,and we know that pi,02 resp. pi-1,02-gi-1,f(euj)(bi-1) are the pressures at the head resp. the foot of the compressor edge euj. So it follows(64)1ujpheuj2=pfeuj2,which is equivalent to the compressor property. This completes the proof.With this theorem we have another characterization of the feasible setM~. The last step in this intersection is the adaption of the results to the spheric-radial decomposition. As mentioned before, the main step is to get a representation of the sets Mi in step (2) of Algorithm 3.Because we assumed that the graph has only one entry, we define the following set:(65)M≔b∈R≥0n∣-1nTb,b∈M^.Consider a sampled point vi∈Sn-1 of step (1) of Algorithm 3 (i=1,…,N). We identify the load vector b with the affine linear function:(66)bir≔rLvi+μ=rωi+μwithωi=Lvi,where μ is the mean value of the Gaussian distribution and L is such that LLT=Σ for positive definite covariance matrix Σ. Note that the index i here is for the sampled vector vi, not for the i-th subgraph. Because bi is defined only for the exit nodes, it must hold b≥0 which leads to the definition of the regular range:(67)Rreg,i≔r≥0∣bir≥0⊇Mi.So the feasible set is(68)Mi=r∈Rreg,i∣birfulfills3,5and7,and with Theorem 5 it holds(69)Mi=r∈Rreg,i∣b~irfulfills25-34,where b~i(r) is from (18). Because b~i depends on r now, the function g depends on r, too. More exactly, g is quadratic in r, so gk(bi(r)) can be represented as follows:(70)r2∑j=1nΨk,j∑l=1nAj,l-1ωi,l2︸a~i+r2∑j=1nΨk,j∑l=1nAj,l-1ωi,l∑l=1nAj,l-1μl︸b~i+∑j=1nΨk,j∑l=1nAj,l-1μl2︸c~i.Here, Ψ is the product of A and Φ. So the inequations, which characterize Mi, are also quadratic in r because the minimal and maximal pressures in the inequalities only count to c~i. Now we can write the feasible sets as follows:(71)Mi=r∈Rreg,i∣0≤ai,tr2+bi,tr+ci,t,where t is an index for the total number of inequalities in Mi. Now the regular range can be cut with the positive intervals of the inequalities and the result is a union of disjoint intervals:(72)Mi=⋃k=1sia_i,k,a¯i,ksi∈N.Here si is the number of disjoint intervals in which all inequalities are fulfilled and a¯i,k and a_i,k are the interval bounds, they depend on the pressure bounds and the controls. So the third step of Algorithm 3 can be executed. With(73)Pb∈M≈1N∑i=1N∑k=1siFχa¯i,k-Fχa_i,k,the probability for a random load vector to be feasible can be computed.


## 4. Existence of Optimal Solutions

In this section we want to have a look at some optimization problems. We distinguish between problems with constant loads (without uncertainty on the demand) and problems with random loads. The latter leads us to optimization problems with so called probabilistic constraints or chance constraints (see [21]). There are two aims in optimization: Minimize the maximal pressure bounds and minimize the controls. For minimizing the maximal pressure bounds, we define the following set:(74)P¯≔⋃i=0npi+,min,∞,in which we look for them. Because the maximal pressures are in proportion to the cost of pipes, an objective function for a optimization problem could be cTp+,max for a cost vector c∈R>0n+1. By minimizing the controls, we use the Euclidean norm of the controls as an objective function.Optimization with Constant Loads.Consider the two optimization problems(75)minp+,max∈P¯cTp+,maxsubjecttob+p+,max∈M^and(76)minu∈R≥1m2u2subjecttob+u∈M^.Lemma 7.
There exists a unique solution for (75).Proof.
We have a look at the side constraintb+(p+,max)∈M^. Using Theorem 5, we can write (75) as(77)minp+,max∈P¯cTp+,maxsubjectto25-35arefulfilled.For i∈{1,…,m2+1} we write inequality (25) (and all remaining analogous) in the following way:(78)0≤gi,kb~i-pi,0min2+pi,kmax2∀k=1,…,niSo every inequality depends only on one upper pressure bound (but more inequalities may depend on the same bound). Because the inequalities are not strict, the values(79)p¯k≔minpkmax∈0,∞∣25-35arefulfilledexist for all k=1,…,n+1 and we can define a value(80)popt,k+,max≔maxp¯k,pkminwhich is obviously a solution of (75). Moreover, this solution is unique, because reducing one component of the upper bounds would either hurt an inequality or it would not fulfill p+,max≥p+,min anymore.For the optimization problem (76) we can formulate a necessary and a sufficient condition for optimal solutions.Lemma 8.
Consider pressure boundsp+,min,p+,max∈Rn+1 s.t. p+,min<p+,max and a suitable vector b+∈Rn+1 (i.e., b~+ fulfills (25)–(27)). For all k=1,…,m2 we set i(k) and j(k) such that f(ek)∈Vi(k)+ and h(ek)∈Vj(k)+.
(a) Ifu∈Rm2 is a solution of (76), it holds for all k=1,…,m2:(81)pik,0min2-gik,fekb~ik≤pjk,0max2.(b) If the inequalities(82)pik,0max2-gik,fekb~ik≥pjk,0min2andpik,0min2-gik,fekb~ik≤pjk,0max2are fulfilled for all k=1,…,m2, an optimal solution u∈R≥1m2 of (76) exists.Note that a vectorb+∈Rn+1 is suitable here if the vector b~+ is feasible for the subgraphs, i.e., the inequalities (25)-(27) are fulfilled for b~+.Proof.
(a) For this part we use a proof by contradiction. Assume that(83)pik,0min2-gik,fekb~ik>pjk,0max2.Because uk≥1 for all k=1,…,m2, it follows(84)pik,0min2-gik,fekb~ik>1ukpjk,0max2.It it Σk∗,i(k)=0 and Πk∗,i(k)=1 because the subgraph with index k∗ is equal to the subgraph with index i(k). And because the subgraphs i(k) and j(k) are neighbours (i.e., directly connected), it follows Πk∗,j(k)=uk and(85)Σk∗,jk=gk∗,jk,1,feuk∗,jk,1bk∗,jk,1=gik,fekb~ik.Thus we have a contradiction to (28) and uk cannot be optimal.
(b) As mentioned before, the vectorb+ is such that b~ is feasible for the subgraphs. So because of the conditions for k∈{1,…,m2} there exists a pi(k),0∈[pi(k),0min,pi(k),0max] and a pj(k)∈[pj(k)min,pj(k)max] so that(86)pik,02-gik,fekb~ik≥pjk,02.With pi(k),02-gi(k),f(ek)(b~i(k))=pf(ek)2 and pj(k),02=ph(ek)2 it follows(87)phekpfek≥1.Therefore we can find a u∈R≥1m2 that the inequalities (28)-(35) are fulfilled. Now we can use theWeierstraß theorem and get the existence of an optimal solution.Note that the sufficient condition in Lemma8 is comparatively strong, there may exist optimal solutions with less strong conditions, e.g., the following.Corollary 9.
If the graph has only one compressor edge, i.e.,EC=1, then it holds: If u∈R≥1 is an optimal solution of (76), u is unique.This conclusion follows directly from the strict monotonicity of the objective function.Optimization with Random Loads.In this section we consider the load vector to be random as we mentioned in Section 2.2. In the optimization problems, we want to make sure that the feasibility of a random load vector is guaranteed for a probability α∈(0,1). This leads us to optimization problems with so called chance or probabilistic constraints. In general problems with probabilistic constraints have the form(88)minfx,ξsubjecttoPhx,ξ≤0≥α,with an objective function f and a vector of uncertainty ξ. There exist many works about probabilistic constraints resp. stochastic programming, e.g., [7] gives an excellent overview about chance constraints in theory and application. In our case we consider the following problem:(89)minp+,max∈P¯cTp+,maxsubjecttoPbp+,max∈M≥α.To handle this problem resp. the probabilistic constraint, the spheric-radial decomposition gives us an explicit representation of the constraint for a Gaussian vector b~N(μ,Σ) with mean value μ and positive definite covariance Σ:(90)Pbp+,max∈M=∫Sn-1μχr≥0∣rLv+μ∈Mdμηv.Further, Algorithm 3 gives a way to approximate this integral in an efficient way. In Section 3 we showed a way to characterize the feasible set M and how to adapt this characterization to the spheric-radial decomposition. This changes our problem with probabilistic constraints (89) to(91)minp+,max∈P¯cTp+,maxsubjectto1N∑i=1N∑t=1stFχa¯i,t-Fχa_i,t≥α.Of course the interval bounds depend on the pressure bounds and the control. For this problem we formulate the next lemma.Lemma 10.
Let a sampling{v1,…,vN} (from step(1) of Algorithm 3) and α∈(0,1) be given. For a suitable μ∈Rn problem (91) has at least one solution p+,max∈P¯.Proof.
For the proof, we fix av∈{v1,…,vN} from the given sampling. From (25) and (26) it follows for all k=1,…,ni with i∈{1,…,m2+1},(92)0≤gi,kb~ir+pi,kmax2-pi,0min2,0≤-gi,kb~ir+pi,0max2-pi,kmin2.With a random vector b~N(μ,Σ)gi,k(b~i(r)) are quadratic functions in r and the upper pressure bounds only influence the constant part with a positive sign:(93)0≤ai,kr2+bi,kr+ci,k+pi,kmax2-pi,0min2,0≤-ai,kr2-bi,kr-ci,k+pi,0max2-pi,kmin2.It follows that a rise in pi,kmax for all k=0,…,ni enlarges the intervals, in which the inequalities hold. The same holds for all inequalities in Theorem 5. The cumulative distribution function of the Chi-Square-Distribution is strictly monotonic increasing and convergent to 1. Thus upper pressure bounds pi,kmax (for i=1,…,n) can be found so that(94)1N∑j=1N∑t=1sjFχa¯j,t-Fχa_j,tis higher that α∈(0,1). The index N is for all sampled points on the sphere and the index sj is for the union of intervals in which all inequalities hold. The interval boundaries a¯j,t and a_j,t depend continuously on the pressure bounds and because the inequalities are not strict, pmax can be found s.t. the side constraint is fulfilled with equality. So there exists at least one solution of (91).


## 5. Numerical Results

In this section we show a few results of implementation. At first we show the idea of the spheric-radial decomposition by using an easy example. Next we show a easy example with one compressor edge and at last we use real data of the Greek gas network. The focus of the implementation is on the theorems proofed in Section4.Example 1.
Our first example is an easy graph without inner control (see Figure6).Assume that node 1 and 2 are gas consumers with mean demandμ=[0.5,0.5]T. The pressure bounds are given by p+,min=[2,1,1]T and p+,max=[3,2,2]T and Φ=E2. Then the incidence matrix is(95)A+=-101-101resp.A=1-101.Consider the inequalities of Theorem 5 to characterize the feasible set. Because the network has no compressor edges the system of inequalities only contains the information (25)-(27), the other inequalities do not occur. The function g is given by(96)gb=10111001·1101b1b2∘1101b1b2=b12+2b1b2+b22b12+2b1b2+2b22.By inserting the pressure bounds, one can see that inequality (25) is always fulfilled and inequality (27) implies(97)b2≤3.With this, inequality (26) follows for 0≤b2≤3:(98)0≤b1≤-b2+8-b22.So the feasible set M is the following (see Figure 7):(99)M=b∈R≥02∣b2≤3∧b1≤-b2+8-b22.Figure 7
Feasible set of Example1.Now consider the optimization problem:(100)minp+,max∈P¯13Tp+,maxsubjecttob∈Mp+,max,with b=μ=[0.5,0.5]T. Lemma 7 guarantees the existence of a unique solution of this problem. We use thefmincon-function inMATLAB®, which minimizes a constrained nonlinear multivariable function. The fmincon default setting uses a interior-point method. For more information we refer to the officialMathWorks® homepage (https://www.mathworks.com/help/optim/ug/fmincon.html). The solver returns the following result:(101)popt+,max=2,3,2.75T.One can easily see that this is the correct solution. The gas demand at node 2 is 0.5, which is equal to the flow q2, so using the equation of momentum implies a difference in the quadratic pressures between node 1 and node 2 of 0.25. The flow through edge 1 is given by the two exit nodes and it is q1=1. This implies a difference in the quadratic pressures between node 0 and node 1 of 1. So the optimal solution popt+,max directly follows.Next we use the spheric-radial decomposition, especially Algorithm3, to calculate the probability P(b∈M) for a random vector b~N(E2,[0.5,0.5]T). Therefore we view the results of 8 tests with 1000 sampled points in each one (Table 1).Table 1
Test results of easy example.
Test 1 Test 2 Test 3 Test 4 Test 5 Test 6 Test 7 Test 8 P(b∈M) 0.3486 0.3481 0.3450 0.3501 0.3484 0.3499 0.3459 0.3475The probability in all eight tests is nearly the same. The mean value is 0.3479 and the variance is2.7723·10-6, which is very small. This shows that the implementation is working nearly exact. The efficiency of the implementation of course is not perfect, but that is not part of this work.Last we want to solve the following optimization problem:(102)minp+,max∈P¯13Tp+,maxsubjectto1N∑i=1N∑t=1sFχa¯i,t-Fχa_i,t≥α.We use again thefmincon-function inMATLAB® to solve the problem eight times, with 1000 sampled points in every case. Because Lemma 8 needs a suitable μ∈R≥02, we assume b~N(0.1·E2,[1.5,1.5]). We set α≔0.9 and the results are shown in Table 2.Table 2
Test results of Example1.
Test 1 Test 2 Test 3 Test 4 Test 5 Test 6 Test 7 Test 8 popt,0max 5.4514 5.4714 5.4758 5.4028 5.4767 5.3921 5.4888 5.5167 popt,1max 2.9802 3.0126 3.0225 2.9725 3.0300 2.9485 3.0336 3.0029 popt,2max 1.0084 1.0004 1.0776 1.0001 1.1221 1.1390 1.0285 1.0558Example 2.
Now we add a compressor edge to the graph of Example1 (see Figure 8).Figure 8
Graph of Example2.Assume that node 1 and 3 are consumers, node 2 is a inner node. The mean demandμ=[0.5,0,0.5]T. The pressure bounds are given by p+,min=[2,1,1,1]T and p+,max=[3,2,2,2]T and ϕe1=1 and ϕe3=1. The incidence matrix is given by(103)A+=-1001-1001-1001resp.A=1-1001-1001.Before we have a look at the inequalities of Theorem 5, we have to remove the compressor edge from the graph. The flows can be computed by using the equation of mass conservation and the new load vector is then(104)b~=b1+q2b2-q2b3withq=b1+b2+b3b2+b3b3.We assume the control to be switched off, so it holdsu=1. With this, the control edge e2 is treated as a flow edge without friction. The resulting graph is shown in Figure 9: from the inequalities (25)–(35) it follows(105)b~3≤3andb~1≤8-b~32,so the feasible set is (see Figure 10):(106)M~=b~+∈Rn+1∣14Tb~+=0∧0≤b~3≤3∧b~2=0∧0≤b~1≤8-b~22.This implies(107)M=b∈R≥03∣0≤b3≤3∧b2=0∧0≤b1≤-b3+8-b32.Figure 9
Separated graph of Example2.Figure 10
Feasible set of the Example2.Solving the problem(108)minp+,max∈P¯14Tp+,maxsubjecttob~∈M~for b=μ=[0.5,0,0.5]T leads us to the optimal solution(109)popt+,max=2,3,3,2.75T.The argumentation for the optimality of this vector is exactly like before in the easy example. Because we have a compressor edge in this graph, we can minimize the control. Therefore we fix the pressure bounds as we did at the beginning of this example and we set b=[1,0,1.5]T. So inequality (31) is not fulfilled anymore for u=1. We use again thefmincon-solver inMATLAB® to solve the following problem:(110)minu≥1u2subjecttob∈Mu.The solver returns the optimal solution uopt=1.1818.Here, one can see that the sufficient condition in Lemma8 is really strong. It’s not fulfilled with the given values but there exists a solution anyway. Of course, the necessary condition is fulfilled.Now, we use again Algorithm3 to compute the probability P(b∈M) for a random load vector b~N(E3,[0.5,0.5,0.5]T). Table 3 shows us the results.Table 3
Test results of control example.
Test 1 Test 2 Test 3 Test 4 Test 5 Test 6 Test 7 Test 8 P(b∈M) 0.1531 0.1500 0.1540 0.1536 0.1530 0.1550 0.1567 0.1528Again we can see, that the probability is nearly the same in all eight cases. The mean probability is0.1535 and the variance is 3.2369·10-6. For the optimization, we reduce the covariance matrix and assume b~N(0.01·E3,[0.5,0.5,0.5]T). For a probability α=0.9 Table 4 states the results.Table 4
Test results of control example.
Test 1 Test 2 Test 3 Test 4 Test 5 Test 6 Test 7 Test 8 popt,0max 3.9013 3.9223 3.9039 3.8901 3.8879 3.9039 3.8901 3.9098 popt,1max 1.6568 1.6565 1.6365 1.6497 1.6496 1.6365 1.6497 1.6492 popt,2max 1.6568 1.6565 1.6365 1.6497 1.6496 1.6365 1.6497 1.6498 popt,3max 1.0000 1.0000 1.0000 1.0000 1.0002 1.0000 1.0000 1.0000Example 3.
In this example, we use a part of a real network. The focus here is on realistic values for our model. We get these values from the GasLib (see [22]), a collection of technical gas network data.“The goal of GasLib is to promote research on gas networks by providing a set of large and realistic benchmark instances.” (http://gaslib.zib.de/)We use a part of the Greek gas network (GasLib-134) (see Figure11).Figure 11
Greek gas network from GasLib [22], red edges are used in our example.For using real values, we first need to know the constantϕe from (5):(111)ϕe=RST2λDLa2.The value RS is the specific gas constant, which can be calculated by RS=R/M with the ideal gas constant R=8.31451kg·m2/s2·mol·K (see [23]) and the molar mass M=16.62kg/kmol (from real data) of the used natural gas. Also the temperature T=289.15K is given from real data. For the (constant) friction factor λ we use the law of Chen (see [24]):(112)1λ=-2log10k/D3.7065-5.0425Relog10k/D1.10982.8257+5.8506Re0.8981.Diameter D=914.4mm and roughness value k=8·10-6m for the first pipe (see Figure 11) are known from real data. The Reynolds-Number Re can be calculated by Re=(vDρ)/η with the dynamic viscosity η=1·10-5Pa·s (see [1]), the mean velocity v=95.48m/s (calculated with the mean volume flow rate which is given by the real data) of the gas and the density ρ=0.7433kg/m3 (from real data). So we get Re=6489590.74≈6.49·106 and with this the pipe friction coefficient after Chen is λ=0.0093.The length of the first pipe is given byL=17.8693km from the real data. Last value missing is the sound speed in natural gas, which we assume to be ideal. The sound speed a depends on the modulus of compressibility (see [9]). For ideal gases, the modulus of compressibility only depends on the adiabatic exponent and the pressure. Using the ideal gas equation simplifies the equation for the sound speed to(113)a=κRTM.The adiabatic exponent for methan, which is the main part of natural gas (>90%) is about 1.33 (see [25]). Together, all values are shown in Table 5 (example for the pipe next to the source in Figure 11).Table 5
Values based on real data.
variable R S T λ D L a ϕ  value 500.27 289.15 0.0093 914.4 17.8693 438.62 1.98 ⋅ 1 0 7  unit m2/s2⋅K K − m m k m m / s m 2 / s 2Now we want to optimize the inner control in the network for a fix load vector. Therefor we choose two nominations:2011-11-01 and2016-02-17. In both cases, thefmincon-function inMATLAB® returns(114)uopt=1as a local minimum, which satisfies the constraint b∈M. This means, for these nominations, a compressor station is not needed so it is switched off. With this, as a last calculation in this paper, we calculate the probability for the given nominations to be feasible. We use the values from the nominations as mean value μ and set b~N(μ·10-3E,μ). The results are shown in Tables 6 and 7.Table 6
Test results for nomination2011-11-01.
Test 1 Test 2 Test 3 Test 4 Test 5 Test 6 Test 7 Test 8 P(b∈M) 0.9849 0.9849 0.9878 0.9885 0.9863 0.9845 0.9854 0.9853Table 7
Test results for nomination2016-02-17.
Test 1 Test 2 Test 3 Test 4 Test 5 Test 6 Test 7 Test 8 P(b∈M) 0.9744 0.9756 0.9788 0.9788 0.9754 0.9775 0.9734 0.9739Again one can see that the results are quite similar. The mean probability in Table6 is 0.9860 and the mean probability in Table 7 is 0.9760. Both values are nearly 1 so there must happen something really unforeseeable that a Greek does not get his demanded gas.


## 6. Conclusion

In this paper, we have used an existing stationary model for gas transport based on the isothermal Euler equations. We have extended this model by compressor stations and presented a characterization of the set of feasible loads. Especially Theorem5 contains the central idea to get an access to the optimization in Section 4 and the spheric-radial decomposition is our main tool to handle the uncertainty in the loads. We were able to show some results for optimization with probabilistic constraints in gas transport. In a next step these results should be generalized.With the characterization of the feasible set, we have shown that the feasibility of a load depends significantly on the pressure bounds. These pressure bounds play a main role in our optimization problems. The other main parts are the inner controls modeled by compressor edges, which should lead in optimal control problems with probabilistic constraints. Here we have built a base for further works with this topic.The idea of separating the graph and removing the compressor edges can be used to extend the model with other components, like, e.g., control valves and resistors. This would change the inequalities in Theorem5, but the approach works also for these modifications. Further, [8] showed a way to deal with a network which is a cycle. Combining both methods, we get an access to model gas networks in high detail.Finally the expansion to a dynamic model might be also possible. The idea is to use a space-mapping approach, because one can use a lot of the analysis we did here. Therefore upcoming papers should deal with these problems and present results of the turnpike-theory, which describes a relation between stationary and instationary solutions. In fact, turnpike-theory implies that the optimal dynamic state is close to the optimal static state in the interior of the time intervals for sufficiently long time horizons. A.J. Zaslavski explains the turnpike-phenomenon in [26]. This book also gives an overview about the theory.



---

*Source: 7984079-2018-09-26.xml*