Emerging Trends on Optimization and Control under Uncertainty in Transportation and ConstructionView this Special Issue
Stationary Gas Networks with Compressor Control and Random Loads: Optimization with Probabilistic Constraints
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.
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  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  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 ), a simplification of the isothermal Euler equation (see e.g., [4, 5]). Further  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 ).
In the next section, we introduce our gas network model, which is based on the model in . 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 Section 3 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 Section 4, 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 Section 5 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 . 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 graph 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 and (). 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 as the set of all flux edges and as the set of compressor edges, it follows with .
Let with denote the balanced load vector and assume for nodes with gas influx and for nodes with gas efflux (). The vector denotes the vector of all ones in the dimension . Furthermore, a vector describes the flows through the edges resp. through the pipes. The pressures at the nodes are defined in a vector . Let pressure bounds be given. For the pressure we consider the constraints . For further modeling we need the following definition:
Definition 1. Consider the graph :(i) denotes the head node of an edge and denotes the feet of an edge for all .(ii) denotes all edges which are connected to node .(iii)The matrix , with is called the incidence matrix of the graph .(iv)For is the (unique) directed path from the root to .(v)For is the (unique) directed path from to .
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
This is equivalent to Kirchhoff’s first law (see ) and by using the definition of the incidence matrix, the equation for mass conservation for the whole graph is
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 byHere, is the gas density, is the speed of sound, is a (constant) friction coefficient and 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 is negligible. These assumptions simplify the equations enormously. With the ideal gas equation (see ) we get a direct dependency between pressure and density. Solving this simplified equation leads us to the so called Weymouth-Equation (see ):Here, is the specific gas constant, is the temperature of the gas and is the length of the pipe . 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].
Compressor stations counteract the pressure drop along the pipes. These stations are modeled as pipes without friction. We use the following model equation (see , section 7.4.7; , section 2.3):Here, and are the pressures at the end resp. at the beginning of one pipe, and are constants and , which is our control, is the specific change in adiabatic enthalpy. More information about the parameters can be found in . In addition,  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:It holds where means that the compressor station is switched off (see Figure 2).
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., . This article also describes, why a multivariate Gaussian distribution is a good choice for random gas demand. So we use some random vector with for a mean vector and a positive definite covariance matrix on a suitable probability space and our aim is to calculate the following probability:for a suitable set .
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 edge gets the number . Our set is defined as follows:
There are many studies for special sets , e.g., for polyhedral sets (see ), 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 called spheric-radial decomposition. It became popular for dealing with probabilistic constraints (e.g., see [8, 14–19]).
Theorem 2 (spheric-radial decomposition). Let be the -dimensional standard Gaussian distribution with zero mean and positive definite correlation matrix . Then, for any Borel measurable subset it holds thatwhere is the -dimensional sphere in , is the uniform distribution on , denotes the -distribution with degrees of freedom and is such that (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 ). Therefore we setand observe that . So our aim is to compute the following probability:So the calculation of the probability (9) for a very general set 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 , Algorithm 4):
Algorithm 3. Let and such that .(1)Sample points uniformly distributed on the sphere .(2)Compute the one-dimensional sets for .(3)Set .
The first step is quite easy. There are many possibilities to sample a set of uniformly distributed points on the sphere , e.g., Monte-Carlo sampling or the randn-generator in MATLAB®. Normalizing these points leads us to the sampling we are looking for. Alternatively, quasi Monte-Carlo methods can be used. In  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 set is convex, these sets are closed intervals but need not to be convex. The aim is to get a representation of () 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 can be calculated easily with the distribution function of the distribution (see  for details of that distribution). With the representation the -distribution can be calculated the following:This leads to the wanted probability. The -distribution function is given by the following (see , p. 132):with , and for all . Note that there are very precise numerical approximations for the distribution function, e.g., the chi2cdf-function in MATLAB® (see Figure 3).
3. Explicit Characterization of the Gas Flow
In this section we want to characterize the feasible set to use Algorithm 3. The aim is to get a representation of the one-dimensional sets of step 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 , Theorem 1. In the proof, the authors use a representation of the pressure loss for the full network graph:which is equivalent to (5) if and only if . In this paper we also consider the case . The matrix is diagonal with the values 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:We want to treat the flows through the compressor edges as loads, so we have to adapt the load vector. The new load vector is defined as follows:The resulting graph is no longer connected (scheme is shown in Figure 4).
We set and with , so consists of subgraphs, which are all connected. Note that is a vector in now. We determine that is the -th self-connected subgraph of (depends on numbering of the graph). Analogously is the -th component of the -th subgraph. The terms and are defined similarly. The feasible set for is then
The last equation in 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 and .
Proof. (i) Consider a vector . It holds and because of (18) it follows directly . After removing the compressor edges, from it follows with (17) that , which is equivalent to . The pressure loss is fulfilled for all flux edges in and because the flow vector doesn’t change in the relevant coefficients while removing the compressor edges, the pressure loss is fulfilled in too. The compressor property is independent of and , so it also holds. From this it follows .
(ii) The other way around is quite similar. We consider a vector . From it follows with (18) that . The equation for mass conservation is fulfilled in every subgraph. So with (17) it follows , which implies the mass conservation in . The pressure loss holds in all self-connected parts of the graph too and because the relevant coefficients of doesn’t change while adding the compressor edges, the pressure loss is fulfilled in too. Last, as said before, the compressor property is independent of and , it is also fulfilled in . This completes the proof.
Because our network graph is a tree, the mass conservation is fulfilled with a full rank incidence matrix of size , so can be computed explicit by the following: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 set . Therefore we will use the following function:which describes the pressure loss between the nodes. We also introduce a new notation:(i) is the -th control on the path from to (ii) is the function for the -th subgraph on the path from to , whereas is the function of the graph . According to this, is the -th component of the function (iii) is the load vector of the -th subgraph in the path from to , whereas is the load vector of the graph and is the -th component of
Analogously we define , and . Furthermore we define the following values:(i) as the largest index of all subgraph, the paths to and pass through(ii) as the number of subgraph which are on the path from from containing and (iii) as the number of controls which are between and
With this notation and these values we define a sum for . If we set and else we setwhere we write only and instead of and for better reading. Last, we define a product andFrom 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.
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 subgraph and node 12 is part of the subgraph , so it follows because node 2 is part of subgraph 2. The notation can be understood as a numbering along the path. It is , the first control on the path from node 2 to node 11 and , the second control on this path. Further it follows () and (). 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 bounds and controls () the following equivalence holds.
A vector with is feasible, i.e., , if and only if the following inequalities hold. For all holds (feasibility inside the subgraphs):For all with holds (feasibility between the subgraphs):
Let us have a look at the inequalities. Inequalities (25)–(27) are well-known from  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 and directly (there may be no direct connection), we have to compare them using the subgraph which is the last part the paths to and pass. Therefore we defined the sums and , which trace back the pressures to .
Remark 6. For an implementation it is wise to define a set containing all controls on the path from to and use an index running over this set to define the sums and the products.
Now we want to proof the theorem.
Proof of Theorem 5. Because the columns of sum up to zero, it holds“⇒” We consider a feasible vector . From the equation of mass conservation of every subgraph it followsWe take the equation of momentum conservation for the -th subgraph () and separate to and :Next we insert (36) into (38) and getWe multiply this equation from left with and insert (37), so it followsand using the definition of the function we haveNote here, that is the vector of pressures of the subgraph and is the vector-valued function, which describes the pressure loss in . Equation (41) holds for all subgraphs with . From this, the inequalities (25)–(27) follow directly with . Next we use (41) component-by-component. For fix with it followsRemember, that is the largest index of all subgraphs, the paths to and pass through. For better reading, we only write instead of . We use the equation of the compressor property for the first compressor on the path from to ,and insert it into (42) for :Note that , and . We use the compressor equation for the second compressor on the path from to and put it, together with (42) into (44):We repeat this procedure for all controls on the path from to . Although we write and instead of and . Thus we getWe take again (42) for the node and do the same procedure just like before for the path from to . With and we getExpanding (46) and (47) and equalizing them (this is possible because ) leads towith resp. defined before. (48) directly implies (28) and (29). Using (48) together with (42) for first time, for second time and for and a third time, the inequalities (30)–(35) follow directly. So this part of the proof is complete.
“⇐” For this part, we consider a vector with , which fulfills the inequalities (25)–(35). We define the following sets ():We setand we will define a value later. So to guarantee such a value exists, we have to show first, that is not empty. Because 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 is not empty. For , Furthermore its because of (25) and (26). The missing intersections are nonempty too, for with , as follows: So the set is not empty. We define the following values:Now we show that our choice of pressures and flows is feasible. It is because of . Further because of for all it holdsandNow we replace in its definition with and because of for all it followsandSo for all and for all it holds which is equivalent toThe equation for mass conservation follows directly from the definition of (). For the pressure loss we useAnd multiply this equation from left with and use (36). With this it follows the equation for momentum conservation:The last equation missing is the compressor property. For a control we fix a vertex () with . We use and solve the formula for in the following way:Consider w.l.o.g. that the subgraphs along the path from to are numbered chronological, s.t. if then . We use to get the same representation for depending on , which is equivalent to . We equalize these representations depending on and and we getWe chose s.t. is the last control on the path from to , so it is . Further we haveand we know that resp. are the pressures at the head resp. the foot of the compressor edge . So it followswhich is equivalent to the compressor property. This completes the proof.
With this theorem we have another characterization of the feasible set . 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 in step of Algorithm 3.
Because we assumed that the graph has only one entry, we define the following set:Consider a sampled point of step of Algorithm 3 (). We identify the load vector with the affine linear function:where is the mean value of the Gaussian distribution and is such that for positive definite covariance matrix . Note that the index here is for the sampled vector , not for the -th subgraph. Because is defined only for the exit nodes, it must hold which leads to the definition of the regular range:So the feasible set isand with Theorem 5 it holdswhere is from (18). Because depends on now, the function depends on , too. More exactly, is quadratic in , so can be represented as follows:Here, is the product of and . So the inequations, which characterize , are also quadratic in because the minimal and maximal pressures in the inequalities only count to . Now we can write the feasible sets as follows:where is an index for the total number of inequalities in . Now the regular range can be cut with the positive intervals of the inequalities and the result is a union of disjoint intervals:Here is the number of disjoint intervals in which all inequalities are fulfilled and and are the interval bounds, they depend on the pressure bounds and the controls. So the third step of Algorithm 3 can be executed. With