Document [original]

Vol.:(0123456789)

Optimization and Engineering

https://doi.org/10.1007/s11081-020-09559-y

1 3

RESEARCH ARTICLE

Controlling transient gas flow inreal‑world pipeline

intersection areas

FelixHennings1,2 · LovisAnderson1 · KaiHoppmann‑Baum1,2 ·

MarkTurner1,2 · ThorstenKoch2,1

Received: 1 November 2019 / Revised: 27 August 2020 / Accepted: 4 September 2020

Abstract

Compressor stations are the heart of every high-pressure gas transport network.

Located at intersection areas of the network, they are contained in huge complex

plants, where they are in combination with valves and regulators responsible for

routing and pushing the gas through the network. Due to their complexity and lack

of data compressor stations are usually dealt with in the scientific literature in a

highly simplified and idealized manner. As part of an ongoing project with one of

Germany’s largest transmission system operators to develop a decision support sys-

tem for their dispatching center, we investigated how to automatize the control of

compressor stations. Each station has to be in a particular configuration, leading in

combination with the other nearby elements to a discrete set of up to 2000 possi-

ble feasible operation modes in the intersection area. Since the desired performance

of the station changes over time, the configuration of the station has to adapt. Our

goal is to minimize the necessary changes in the overall operation modes and related

elements over time while fulfilling a preset performance envelope or demand sce-

nario. This article describes the chosen model and the implemented mixed-integer

programming based algorithms to tackle this challenge. By presenting extensive

computational results on real-world data, we demonstrate the performance of our

approach.

Keywords Transient gas network optimization· Mixed-integer programming·

Realistic modeling of compressor stations

* Felix Hennings

[email protected]

Extended author information available on the last page of the article

F.Hennings et al.

1 3

1 Introduction

Throughout the past years, the mathematics of gas transport has been an intensively

studied topic. Natural gas was, is, and will be one of the primary energy sources in

Germany, making efficient and safe transport a field of high economic and politi-

cal relevance according to the German Federal Ministry for Economic Affairs and

Energy (2019). Furthermore, the task is also challenging from a mathematical point

of view. One such challenge is in the compressors, which push the gas through the

network by increasing its pressure. Compressors are typically set up as a compres-

sor station consisting of multiple compressor units, each representing the union of

one single compressor machine and its associated drive. These compressor units are

dynamically switched on or off and used in different sequences to meet the current

needs in terms of compression ratios and flow rates. Out of the theoretically pos-

sible arrangements of compressor units, the set of technically feasible arrangements

are known as the configurations of a compressor station. The intersections of major

transportation pipelines are usually composed of multiple compressor stations as

well as other elements like valves or pressure regulators. Such a structure makes it

possible to adjust the connections of the intersecting pipelines and operate the sys-

tem for various pressure levels. To optimize control of these areas, which includes

all the discrete decisions that need to be made for the single elements while taking

into account the corresponding technical restrictions, is itself already a combinato-

rial challenge. The physics of gas flow further increases the complexity of the prob-

lem. In pipes, the Euler Equations describe the flow of gas. This set of nonlinear

hyperbolic partial differential equations (PDEs) yields even in simplified versions

computationally challenging constraints, see Osiadacz (1996).

Research in this field first focused on the simulation of gas flow, i.e., dealing

with the partial differential equations given all the discrete decisions. This topic

has been studied for many decades already, see for example Brouwer etal. (2011)

and the references therein. In recent years however, the optimization of gas

transport, including the consideration of the combinatorial aspects, has gained

increased attention. In the article of Ríos-Mercado and Borraz-Sánchez (2015), a

general overview of optimization problems related to natural gas is given, includ-

ing the transport of gas. The majority of the referenced literature considered

the stationary (or steady-state) gas transport problem. There, they deal with one

stable network state which allows for an algebraic description of the gas flow.

Koch etal. (2015) and Pfetsch etal. (2014) give an overview of state-of-the-art

approaches tackling the stationary problem on large real-world instances with a

considerable amount of detail regarding the different network elements.

The more challenging variant is the transient gas transportation problem, which

we consider in this article. Here, the goal is to find a set of control decisions for the

elements over a future time horizon discretized into individual time steps. For this

problem, research is still in the early stages. One of the first publications on transient

gas transport optimization was the thesis of Moritz (2007) in which they presented

a mixed-integer programming (MIP) model for the problem. In contrast to the pos-

sible configurations described above, they only used a model consisting of single

1 3

Controlling transient gas flow inreal‑world pipeline…

compressor units, whose compression capabilities are limited by a given minimum

and maximum value for the power used for the compression. The nonlinear con-

straints resulting from these power bounds, as well as the nonlinear pipe equations,

are approximated by piecewise linear functions. To solve the model for the objective

of minimizing a combination of the compressors’ fuel consumption and the cost to

switch the elements, they used a special branching scheme for the piecewise lin-

ear functions as well as a simulated annealing heuristic described in the article of

Mahlke etal. (2007). A little later, Domschke etal. (2011) also presented a solution

to the problem of minimizing the fuel consumption of compressors. They consid-

ered a similar problem formulation as Moritz (2007), using the same model for com-

pressors and approximating the nonlinear constraints by piecewise linear functions.

However, they combined solving this MIP with solving a nonlinear problem (NLP)

formulation of the problem in an alternating way. From the solution of that nonlin-

ear problem, they deduce a refinement of the piecewise linear approximations and

repeat this procedure until finally ending up with a solution to the overall mixed-

integer nonlinear problem (MINLP) within a chosen approximation error. In other

articles dealing with transient transport optimization, the discrete nature of the prob-

lem is neglected, which leads to NLP formulations. As an example, we mention the

work of Zlotnik etal. (2015) and Mak etal. (2016), who decide on the compression

ratios of compressors that are assumed to be always active, while again minimizing

their fuel consumption. Very recently, a few more studies on transient gas network

optimization have been published. Hahn etal. (2017) use a specialized branching

rule to solve a MINLP formulation of the problem with the objective of minimizing

fuel consumption. For the compressors, they introduced the theoretical concept of

different modes to switch between configurations, which each have separate feasible

regions. However, in the end, they restrict to exactly one mode with nearly unlim-

ited compression capabilities for their experiments. Another approach combining

different specialized solving techniques is presented by Gugat etal. (2018), where

a MIP model and an NLP model are iteratively solved for each single time step to

determine an instantaneous control of the problem, i.e., a control that only looks

one time step ahead when making a decision. The two used models are based on

a specific discretization of the Euler equations. For compression, both models use

single compressor units featuring a linear feasible region. In contrast to the other

previously mentioned publications, Gugat etal. (2018) do not minimize fuel con-

sumption. Rather, the objective is to deviate as little as possible from a set of future

pressure and flow values given at the boundary nodes. To close, we mention the

work of Burlacu etal. (2019), who consider maximizing the amount of temporarily

stored gas in the network while maintaining a feasible transient control of the ele-

ments. They also introduce a new discretization of the Euler Equations resulting in

a formulation close to the algebraic form of the stationary model. However, for the

compressors, they model only single units restricted by upper and lower bounds in

the compression ratio as well as the absolute pressure difference. This problem is

solved by alternating between solving a MIP model, which they obtain by replacing

the nonlinear constraints by piecewise linear functions, and an NLP model, in which

the discrete decisions from the MIP are already fixed. Using this approach they are

able to obtain globally optimal solutions for their model.

F.Hennings et al.

1 3

All approaches in the publications on transient gas network optimization men-

tioned above use a somewhat idealized model for the compressor stations. In con-

trast to this, we present in this article a transient gas network optimization problem

featuring a compressor model with a so-far unmatched amount of detail. The model

is based on the modeling presented by Koch etal. (2015) for the stationary case and

takes the above described substructures for compressors into account. We consider

the problem not on the complete network, but on those individual network areas that

contain the compressor stations as well as additional active elements. We call these

subgraphs network stations, an example station is presented in Fig.1. The set of

all network stations contains the majority of active elements in the whole network.

Regarding the number of contained elements, a single network station is comparable

to or even larger than the networks considered in the literature on transient optimi-

zation. Hence, they constitute challenging instances even though only representing

a subgraph of the original gas network. Furthermore, precise control recommenda-

tions for network stations can help gas network operators to manage critical situa-

tions that arise during day-to-day operation.

One difference to general gas network problems is the shortness of pipes due to

the proximity of the elements in a network station, see Table1 in Sect.5. Because of

their shortness, the pipes’ ability to store gas is negligible owing to their small vol-

umes. Furthermore, the pressure loss induced by friction in the pipe depends on its

length and has therefore reduced impact. This allows us to use a linear pipe model

as introduced by Hennings (2018) without losing much accuracy and still producing

realistic results.

As inputs to the problem, we are given the topology of the network station under

consideration, an initial state for each of the contained network elements, and future

demands in terms of both inflow and pressure levels at the station’s boundaries. The

goal is to then find a feasible control of all the network elements over time, which

can be interpreted as a recommendation for network operators on how to best oper-

ate the network station in the future. The overall objective is split into two parts and

will be realized as the minimization of one weighted sum over the different terms.

On the one hand, we try to meet the future demands as best as possible, i.e., allow-

ing and penalizing the deviation from these. On the other hand, we minimize the

total number of control changes. Having such a stable control with a small number

of changes is preferable since it reduces strain on the technical elements and enables

the gas network operators to understand and actually perform the desired control

recommendations.

To solve this problem, we formulate it as a MIP model. However, the model

turns out to be quite challenging, according to our experiments. On larger net-

work stations, state-of-the-art MIP solvers cannot solve it reliably, even for a

small number of future time steps. For this reason, we present a heuristic algo-

rithm that can efficiently and reliably solve the problem also on harder instances.

Its main idea is to determine a majority of the discrete control decisions by solv-

ing independent stationary variants of the original MIP model for each of the

future time steps. Based on these we create a solution for the overall transient

problem by again solving a variant of the original MIP. Apart from the complex-

ity reduction gained from splitting the MIP model, our algorithm also allows to

1 3

Controlling transient gas flow inreal‑world pipeline…

enforce certain technical restrictions outside of its MIP subroutines. Note that the

decisions determined by the stationary model yield good solutions for the origi-

nal transient problem too, since network stations only have a limited amount of

pipeline volume due to the pipes’ shortness. Therefore, only small amounts of

gas can be stored in the network station, which makes the corresponding transient

optimization problem a “stationary-flavored” one.

The rest of the article is structured as follows: In Sect.2, we will describe the

mathematical constraints for all used elements and formulate the original MIP

model. We highlight that we do not present a MIP formulation for the restrictions

introduced in Sect.2.8 since these will be satisfied by our heuristic algorithm with-

out including them in the solved model variants. In Sect.3, we describe the preproc-

essing needed to convert the given compressor data into a linear description of the

feasible operating range. Our heuristic algorithm is the topic of Sect.4, which also

introduces additional variants of the MIP model. Sect.5 then presents the results of

our algorithm when applied to a large number of real-world network situations. We

finish by concluding this work in Sect.6.

2 Mathematical model

In this section, we give a detailed description of the considered transient gas network

optimization problem on a network station. At the same time, we introduce the cor-

responding variables and constraints to create a MIP formulation for the problem.

Fig. 1 Example of the medium-size network station E, see Table1 for more details to its properties. The

colored triangles represent the entry and exit nodes of the station. Furthermore, we have denoted

the single network elements by (pipe), (valve), (regulator), and (compressor station).

(Color figure online)

F.Hennings et al.

1 3

The exception to this is the restrictions presented in Sect.2.8 for which the formu-

lation is not given. Our algorithm developed in Sect.4 produces control decisions

that satisfy these restrictions without their explicit inclusion in any MIP subroutine.

Therefore, a corresponding MIP formulation is not required for these restrictions.

We model the topology of a network station as a directed graph

G=(V,A)

which the arcs

represent the different network elements and nodes

represent the

junctions of the arcs. We split

into individual sets

A=Api ∪Ava ∪Ars ∪Arg ∪Acs

for the network elements contained in a network station, i.e., pipes, valves, resistors,

regulators, and compressor stations. Note that regulators are also often referred to

as control valves in literature, e.g., in Ehrhardt and Steinbach (2003) or Koch etal.

(2015). Likewise, we split the set of nodes

V=Vb∪V0

into boundary nodes and

inner nodes of the network station, respectively. Here, boundary nodes

repre-

sent those having inflow and pressure level demand values for future time steps. We

define the set of considered time steps as

T0∶= {0, …,k}

where

T∶= T0⧵{0}

the set of future time steps having demand conditions at boundary nodes. Each time

step t has an associated value

𝜏(t)

representing the time difference in seconds from

t to the initial state time step 0. However, some restrictions may not be aligned with

this discrete set of time points but consider intermediate time points instead. For

those cases, we define that all variable values at some

t∈T0

are present for the open

interval between t and the subsequent discrete time point

t+1

at which a new value

might be set.

The essential quantities we use to describe the gas flow are pressure and mass

flow. We have pressure variables

pv,t

at each node

v∈V

and time

t∈T0

as well as

variables

qa,t

for the flow from

𝓁

to r on each non-pipe arc

(

𝓁

,r)=a∈

⧵

and

time

t∈T0

. For a pipe

(

𝓁

,r)=a∈

, we have two flow variables

q𝓁,a,t

represent-

ing the inflow into the pipe at end node

𝓁

and

qr,a,t

representing the outflow out

of the pipe at end node r, similar to the model of Ehrhardt and Steinbach (2003)

on pipes or the one of Domschke etal. (2011) on all arcs. The last flow quantity

we introduce is the inflow

dv,t

entering the network station at each boundary node

v∈Vb

and

t∈T0

. Although we are given flow demands for the future, we allow

deviations from them, and hence have to have a variable capturing the actual inflow

value.

We assume that there are upper and lower bounds on all stated quantities for each

point in time

t∈T0

, represented by bars above and below the variable, respectively.

Note that while pressure is always positive, flow can be negative, representing flow

against its orientation in the case of arcs and flow out of the network in the case of

boundary nodes.

Gas flowing through a network always induces a loss of pressure in flow direction

due to friction. However, we only consider friction-induced pressure loss in pipes

and resistors. For valves, regulators, and compressor stations, we either assume the

friction to be small enough to be negligible or to be captured by resistors surround-

ing the element in the network station.

For the formulation of the MIP model, we often use indicator constraints. Here,

a constraint

x≤a

for a set of variables

x={x1,…,xi}

with

a∈ℝi

a0∈ℝ

only active when a logical condition holds, i.e., a binary variable y attains the value

b∈{0, 1}

. We denote such a constraint by

1 3

Controlling transient gas flow inreal‑world pipeline…

Since we have bounds on all the variables, constraints of this kind can be reformu-

lated as linear constraints in a bigM approach, see Bonami etal. (2015) for more

details. A list of all used variables can be found in the appendix in Table3.

2.1 Compressor stations

Compressor stations are the most complex elements in the network. They are

responsible for increasing the pressure and thereby are essential in controlling the

gas flow. We based our model on the description of Fügenschuh etal. (2015), where

the elements are called compressor groups instead.

2.1.1 Structure ofacompressor station

A compressor station

(𝓁,r)=a∈Acs

has three different modes, referred to as

bypass, closed, and active mode. In bypass mode, the element is bypassed and there-

fore allows unrestricted gas flow without changing the pressure level. For the closed

mode, the element blocks the gas flow, which disconnects the network between both

end nodes. Finally, in active mode, gas is compressed and the pressure is increased

along the direction of the flow.

When compressing, the station uses its set of associated compressor units

. These represent the union of a single turbo compressor machine, the tech-

nical element performing the actual compression of the gas, and a correspond-

ing drive delivering the power needed for this process. In the compressor station,

these units are combined in series and/or parallel to allow proper reactions to dif-

ferent compression requirements. The set of all allowed serial-parallel compres-

sor unit combinations is called the set of configurations

for a compressor sta-

tion

a∈Acs

, from which exactly one has to be chosen if the compressor station

is in active mode. For each of these configurations

c∈Ca

, we create a polytope

describing the feasible operating range of the compressor station using configura-

tion c. The polytope is embedded in the three dimensional space of incoming pres-

sure

p𝓁

, outgoing pressure

and flow

. It is described as the intersection of a

set of half-spaces

= {(𝛼

,𝛼

)∈ℝ

encoding inequalities of the form

𝛼0⋅p𝓁+𝛼1⋅pr+𝛼2⋅q+𝛼3

≤

. In Sect.3 we describe how we create these fea-

sible operating range polytopes for the single configurations of the compressor sta-

tions based on the operating ranges of the corresponding compressor units used in

the configuration.

2.1.2 Compressor station model

To model the above-described constraints, we use the disjunctive formulation of

Balas (1985, 2018). Thus, we introduce binary “selection” variables

mcf

c,a,t

for each

configuration

c∈Ca

of a compressor station

(𝓁,r)=a∈Acs

, as well as

mby

a,t

and

mcl

a,t

for bypass and closed mode, respectively, where a value of 1 indicates the selected

b→a

x≤a

F.Hennings et al.

1 3

configuration or mode. Furthermore, we introduce corresponding sets of pressure

and flow variables:

Note that we do not need all variables for each mode since

p𝓁=pr

holds true for

bypass, and

q=0

holds true for closed. Using the binary variables, we force the

pressure and flow variables of all non-selected configurations or modes to be zero

and only enforce the constraints of the selected configuration or mode. The actual

pressure and flow variables of the compressor station can hence be represented

as the sum of the corresponding variables for the configurations and modes. The

bounds of these variables specific to the configurations and modes are equal to those

for the compressor station. The bound interval may not include zero however. In

these cases we expand the interval by replacing its endpoint with the smallest abso-

lute value by zero to make the non-selected variable values feasible.

The formal constraints for all

(𝓁,r)=a∈Acs

and

t∈T

are then given as:

pby

a,tq

a,t

bypass mode variables

l-cl

a,tpr-cl

a,t

closed mode variables

l-cf

c,a,t

pr-cf

c,a,t

qcf

c,a,t

∀c∈C

aconfiguration variables

(1)

∑

c∈

mcf

c,a,t+m

a,t+mcl

a,t

(2)

p𝓁,t=p

a,t+pl-cl

a,t+

∑

c∈C

pl-cf

c,a,t

(3)

pr,t=p

a,t+pr-cl

a,t+

∑

c∈Ca

pr-cf

c,a,

(4)

a,t=q

a,t+

∑

c∈

qcf

c,a,

(5)

l-cf

c,a,t

c,a,t≤p

l-cf

c,a,t≤p

l-cf

c,a,tm

c,a,t∀c∈C

(6)

r-cf

c,a,t

c,a,t≤p

r-cf

c,a,t≤p

r-cf

c,a,tm

c,a,t∀c∈C

(7)

qcf

c,a,t

c,a,t≤q

c,a,t≤q

c,a,tm

c,a,t∀c∈C

(8)

pby

a,t

a,t≤p

a,t≤p

a,tm

1 3

Controlling transient gas flow inreal‑world pipeline…

2.2 Pipes

Gas flow in pipelines is for operational purposes modeled as one-dimensional flow

through a straight cylindrical pipe. We use the friction dominated model derived from

the isothermal Euler Equations and the equation of state for real gases, see the model

(FD1) of Brouwer etal. (2011) and the model (ISO3) of Domschke etal. (2017).

For the discretization, we use the implicit box scheme introduced by Domschke etal.

(2011); Kolb etal. (2010). Here, the spatial domain is the length

of pipe

a=(𝓁,r)

and the time domain is the set of time steps

. For two adjacent time points

and

the discretized model reads as

Here

denotes the specific gas constant,

the gas temperature,

the compressibil-

ity factor,

the cross-sectional area of the pipe,

𝜆a

the friction factor of the pipe,

the diameter of the pipe, and

the gravitational acceleration. With

sa∈[−1, 1]

denote the slope of the pipe, i.e., the quotient of the elevation increase between the

pipes endpoints and the length

of the pipe. We assume that all these values are

predetermined parameters.

For the friction factor, we use the formula of Nikuradse (1950) for hydraulic

rough pipelines assuming full turbulence:

(9)

a,t

a,t≤q

a,t≤q

a,tm

(10)

l-cl

a,t

a,t≤p

l-cl

a,t≤p

l-cl

a,tm

(11)

r-cl

a,t

a,t≤p

r-cl

a,t≤p

r-cl

a,tm

(12)

𝛼

l-cf

c,a,t+𝛼1p

r-cf

c,a,t+𝛼2q

c,a,t+𝛼3m

c,a,t≤0

∀(𝛼0,𝛼1,𝛼2,𝛼3)∈Hc∀c∈C

mby

a,t,mcl

a,t∈{0, 1}

mcf

c,a,t

∈{0, 1}∀c∈C

(13)

p𝓁,t1+pr,t1−p𝓁,t0−pr,t0+

(𝜏(t

)−𝜏(t

))

LaAa

(

qr,a,t1−q𝓁,a,t1

)

(14)

pr,t1

−p𝓁,t1

𝜆aRsTzaLa

4DaA2

(

|q𝓁,a,t1|q𝓁,a,t1

p𝓁,t1

+|qr,a,t1|qr,a,t1

pr,t1

)

+gsaLa

p𝓁,t1

+pr,t1

)

F.Hennings et al.

1 3

Apart from the pipe diameter, it depends only on the integral roughness parameter

of the pipe and is commonly used in gas network optimization. As an example,

we refer to the article of Pfetsch etal. (2014), which furthermore provides additional

background information on friction factor formulas. Regarding the compressibil-

ity factor, the approximation formula developed by Papay (1968) is used, see Saleh

(2002) for an English reference. It is valid for natural gas of up to 150bar and there-

fore fits our use-case of pressure values usually being between 1bar and 100bar.

The formula is based on the pseudo-critical pressure

and temperature

. Both

values depend on the mixture of the gas, which we assume to be constant all over

the network. The constant compressibility factor

per pipe

a=(𝓁,r)

, which we

use in the pipe equations, is determined as the average of the corresponding val-

ues derived from the initial state pressure values of the pipes’ end nodes, i.e., by

za=(

(

p𝓁

,0)+

(

r,0))∕2

In our model, we use a linearized variant of Eq.(14) as proposed by Hennings

(2018). Analogous to them, we fix the absolute velocity

|q|

to a predefined

constant in the friction term. As mentioned before, we found that pipes in network

stations are relatively short in comparison to those in the rest of the network. For

example, the network stations considered in our computational experiments in

Sect.5 have an average length per station of at most 500m, according to Table1. On

the other hand, the transport pipelines that are incident to the stations are routinely

100km or longer. Therefore, the pipes have much less impact inside stations, which

allows us to assume a linear pipe model without losing much in terms of overall

accuracy.

The final equations for each pipe

(

𝓁

,r)=a∈Api

and all adjacent time points

t0,t1∈T0

reduce to:

Note that (13) stays the same as above. The constant velocity

|vx,a|

for one of the end

nodes

x∈{𝓁,r}

of pipe

a=(𝓁,r)

is determined based on the flow and pressure val-

ues of the given initial state, i.e.,

𝜆

(

2 log10

(

+1.138

)−2

(p)=1−3.52 p

pce−2.26 T

Tc+0.247

(

e−1.878 T

𝓁,t1+pr,t1−p𝓁,t0−pr,t0+

(𝜏(t

)−𝜏(t

))

LaAa

(

qr,a,t1−q𝓁,a,t1

)

=0(13

)

(15)

pr,t1

−p𝓁,t1

𝜆

4DaAa

(

|v𝓁,a|q𝓁,a,t1

+|vr,a|qr,a,t1

)

+gsaLa

(

p𝓁,t1

+pr,t1

)

1 3

Controlling transient gas flow inreal‑world pipeline…

2.3 Resistors

Resistors are artificial elements created to model places of high friction in the net-

work, e.g., caused by non-standard network elements without further relevance for

the operation of the network. The pressure drop over a resistor arc

(

𝓁

,r)=a∈Ars

at time

t∈T

is, according to Fügenschuh etal. (2015), defined by the Darcy-Weis-

bach equation:

Here, the friction factor of the resistor is called the drag factor and is represented

by the parameter

𝜁a

. We determine the constant compressibility factor

in the same

way as for pipes. Note that the formula depends on the flow’s direction due to

pin,t

representing the pressure of the incoming gas. The similar looking friction term in

the Eqs.(14) and (15) of the pipe model does not contain this dependency thanks to

the chosen discretization scheme.

Analogous to the approach for pipes in Eq. (15), we linearize the model by

assuming a constant velocity

|q|

, which also includes the flow direction

dependent pressure value. The equation for each arc

(𝓁,r)=a∈Ars

and time

t∈T

then reads as

The constant velocity value is again calculated based on the initial element state

and is defined as an average of the two velocities using the pressure from the cor-

responding resistor’s end nodes as

2.4 Valves

By being open or closed, valves connect or disconnect two nodes of the network.

Closed valves work exactly as described by the closed mode of a compressor sta-

tion, while open valves imply the same behavior as the corresponding compres-

sor station’s bypass mode. The valve’s mode is captured by the binary variable

mop

a,t

, which attains the value 1 for the open mode and 0 for closed. For a valve arc

(

𝓁

,r)=a∈Ava

and time

t∈T

, we can capture the modes’ behavior using indicator

constraints as:

vx,a

RsTza

x,a,0

x,0

𝓁,t−pr,t=

𝜁aRsTza

2A2

a(|q

a,t

in,t).

(16)

p𝓁,t−pr,t=

𝜁

qa,t

∶= 1

(

RsTza

a|q

a,0

p𝓁

+RsTza

a|q

a,0

r,0 ).

F.Hennings et al.

1 3

2.5 Regulators

A regulator or control valve is a valve with a variable opening degree, used to

reduce the pressure along the direction of the flow. Regulators can reduce the pres-

sure in active mode and also be bypassed or closed like a compressor station. For

each mode, there is a binary selection variable, where a value of 1 again indicates

the uniquely selected mode per time. For each arc

a∈Arg

and all times

t∈T

, the

behavior is captured by the following constraints:

Note that for the logical condition in (23), the sum

mby

a,t+mac

a,t

can only be

∈{0, 1}

due to (21) and therefore acts as an implicit binary variable. Furthermore, all reg-

ulators have a flap trap preventing gas flowing against the topological orientation,

which is modeled by (25).

2.6 Nodes

The nodes do not represent technical elements but rather establish the connec-

tions between them. The pressure coupling is realized by using the node’s pressure

(17)

mop

a,t

=1→p𝓁

≤p

r,t

(18)

mop

a,t

=1→p𝓁

≥p

r,t

(19)

mop

a,t

=0→q

a,t

≤

(20)

mop

a,t=0→qa,t≥

mop

a,t

∈{0, 1}.

(21)

a,t

+mac

a,t

(22)

mby

a,t

=1→p𝓁

≤p

r,t

(23)

mby

a,t

+mac

a,t

=1→p𝓁,t≥pr,

(24)

mcl

a,t

=1→q

≤

(25)

≤

a,t

,mby

a,t

,mac

a,t

∈{0, 1}

1 3

Controlling transient gas flow inreal‑world pipeline…

variables in the constraints of all its incident arcs. To connect the mass flow vari-

ables of the arcs, we have the following flow conservation constraints for all

t∈T

Note that pipes have separated inflow and outflow variables, while all other arcs

only have a single flow variable, as explained in the beginning of Sect.2.

2.7 Network station

In addition to the constraints imposed by the single elements contained in the net-

work station, we also consider restrictions associated with the network station itself.

The model described below is an extension of the one described in Fügenschuh etal.

(2015) for similar structures. There, network stations are called compressor stations,

and the elements we call compressor stations are called compressor groups. Our

project partner suggested adjusting the naming scheme here since the previous one

can be misleading, especially so as a network station does not need to contain an

element with compression capabilities.

2.7.1 Network station structure

For a network station, two main decisions have to be made for each point in time.

We have to choose exactly one flow direction from the set

of available flow direc-

tions and exactly one operation mode from the set

of available operation modes.

A flow direction represents the general orientation of the gas flow through the

station, for example, flow from the northern and eastern boundary nodes to those at

the southern border of the station. Hence, a flow direction is defined as a partition

of the boundary nodes into entry nodes, exits nodes, and no-flow nodes, at which

respectively, gas enters the station, gas leaves the station, or the gas inflow and out-

flow is zero. Only a subset of the theoretically possible partitions is feasible for the

station, forming the set

Similarly, the modes and configurations of the network station’s elements cannot

be arbitrarily combined. For example, separating an active compressor station from

the network by closing all surrounding valves would not be feasible in reality. There-

fore we are given the operation modes, which represent the available mode and con-

figuration combinations of all valves and compressor stations of the network station.

(26)

∑

(𝓁,v)=a∈Api

qv,a,t−

∑

(v,r)=a∈Api

qv,a,t

+∑

(

𝓁

,v)=a∈A⧵A

qa,t−

∑

(v,r)=a∈A⧵A

qa,t+dv,t=0∀v∈V

(27)

∑

(𝓁,v)=a∈Api

qv,a,t−

∑

(v,r)=a∈Api

qv,a,t

+∑

(

𝓁

,v)=a∈

⧵

Api

qa,t−

∑

(v,r)=a∈

⧵

Api

qa,t=0∀v∈V0

F.Hennings et al.

1 3

The set of feasible mode combinations for valves and compressor stations is, most

of the time, very small in comparison to the set of all possibilities. Often, a feasible

combination prescribes only one or a few specific paths for the gas flow through the

station using non-closed elements. However, imagine a situation in which a com-

pressor station is in active mode in a particular operation mode of the network sta-

tion. In this case, all available configurations of the compressor station generally

yield feasible operation modes when combined with the same settings for the other

elements.

Note that operation modes do not prescribe the mode of regulators. For these, the

network operator cannot set the desired mode directly but chooses target values for

the incoming and outgoing pressure of the element as well as its flow. The regula-

tor tries to meet these values by adjusting its opening degree according to a given

set of rules, which also determines the mode. Since target value modeling is out

of scope for this paper, we simply do not include regulator modes in the operation

mode restrictions.

Finally, not all combinations of flow directions and operation modes are valid.

One could, for example, think of a compressor station that is only active when the

network station operates in a south to north flow direction. Hence, we are given the

set

Q⊂F×O

of valid combinations of flow directions and operation modes.

2.7.2 Operation modes model

We introduce binary variables

mom

o,t

for each

o∈O

and

t∈T0

, where a value of 1

represents the selection of o at time t. Furthermore, we define the function M(o,a)

that maps operation modes to the individual modes and configurations of valves and

compressor stations:

Using M(o,a), we can state the constraints modeling the relations related to opera-

tion modes for all

t∈T

as:

M(o,a) ∶= xwhere xis the mode or configuration of arc a

in operation mode o∀o∈O∀a∈Ava ∪A

with x∈{op, cl}if a∈Ava

x∈{by, cl}∪

aif a∈

Acs

(28)

∑

o∈O

mom

o,t=

(29)

mop

a,t=

∑

o∈

∶M(o,a)=op

mom

o,t∀a∈A

(30)

mby

a,t=

∑

o∈O∶M(o,a)=by

mom

o,t∀a∈A

1 3

Controlling transient gas flow inreal‑world pipeline…

Note that either Eqs.(30), (31) or one of the constraints of Eq.(32) for one of the

configurations

c∈Ca

can be omitted, since it follows from the remaining constraints

combined with Eq.(1).

2.7.3 Operation mode unavailability

Certain operation modes are not available at specific points in time, and we set the

corresponding variables

mom

o,t

for each associated operation mode o and time t to

zero. The basis for this is the non-availability of compressor units over time, which

is part of the model input data. As explained in Sect.2.1, a configuration

c∈Ca

some compressor station

a∈Acs

represents the serial and/or parallel combination of

a subset of the compressor station’s compressor units. Hence, the unavailability of

a certain compressor unit at time t results in the unavailability of all configurations

using this unit, which causes all operation modes using these configurations to be

unavailable as well.

The unavailability period of a compressor unit may not be aligned with the set

of discrete time points

. According to the definition given at the start of Sect.2,

the network station operation mode of some time

t∈T0

is also used in the open

interval

(t,t+1)

. Thus, the unavailability of an operation mode

k∈(t,t+1)

with

t,t+1∈T0

results in the unavailability of that operation mode for time t. At the

same time, it stays potentially available for time

t+1

2.7.4 Flow directions model

Similar to the way we handle operations modes, we introduce binary variables

mfd

representing the selection of a flow direction

f∈F

at time

t∈T0

if the variable

attains the value 1. As explained above, a flow direction

is defined by the signs of

the flows on boundary nodes. We write

for entry nodes with positive inflow and

f−

for exits with negative inflow, i.e. outflow. Hence, using the power set

, a flow

direction

is defined as

Note that if

v∉f+

and

v∉f−

for flow direction

(f+

,f−)

, then node v is a no-flow

node, and its inflow and outflow are zero.

Using the inflow variable

dv,t

for boundary node v and time t, we can define the

flow direction constraints for each

t∈T

as:

(31)

a,t=

∑

o∈O∶M(o,a)=cl

mom

o,t∀a∈A

(32)

c,a,t=

∑

o∈O∶M(o,a)=c

mom

o,t∀c∈Ca∀a∈A

mom

o,t

∈{0, 1}∀o∈O.

(f+,f−)=f∈

⊆

(

b)×

(

b)where f+∩f−=�.

F.Hennings et al.

1 3

Note that in the constraints(34) and (35), the sum of the

mfd

f,t

variables act as an

implicit binary variable since its value can only be in

{0, 1}

due to (33).

The restriction to feasible combinations of operation modes and flow directions is

implemented as:

2.7.5 Flow direction exit pressures

Apart from the consequences that the flow direction choice has on the correspond-

ing boundary node inflows, it may add restrictions to the upper pressure bounds on

some boundary nodes. Each node v in the set

Vb-ex ⊆Vb

is given a limiting upper

pressure value p

exit

that is only considered if the node is in the outflow set of the

currently active flow direction, i.e., serving as exit node of this flow direction. These

exit pressure restrictions originate from the historical development of the network

infrastructure and may, for example, occur when a transport pipeline has been

upgraded to allow for higher maximum pressure values. If then the adjacent network

station contains equipment, for example a metering station, that is a) only used when

gas leaves the station in the direction of this pipeline and b) has an upper pressure

bound in line with the old pipeline, we have to model this bound explicitly when gas

is flowing in this direction. Before the pipeline’s upgrade, this upper pressure bound

was implied by the corresponding pressure bounds of the nodes in the pipeline. The

corresponding constraint is the following for each time

t∈T

Again, the sum of the

mfd

f,t

variables acts as an implicit binary variable due to (33).

2.7.6 Flow direction conditions

Denoting the final constraints concerning flow directions, there can exist a special

set of conditions

that concern the amount of flow over sets of boundary nodes.

Each condition

w=(f,Vw

,Vw

)∈W

states that the flow over the set of boundary

(33)

∑

f∈F

mfd

f,t=

(34)

∑

f=(f

−

)∈

∶v∉f

−

mfd

f,t=1→0≤dv,t∀v∈V

(35)

∑

f=(f+,f−)∈F∶v∉f+

mfd

f,t=1→0≥dv,t∀v∈V

mfd

f,t

∈{0, 1}∀f∈F.

(36)

o,t≤

∑

(f,o)∈

mfd

f,t∀o∈Ot∈T

(37)

∑

f=(f

−

)∈F∶v∈f

−

mfd

f,t=1→pv,t≤pexit

v∀v∈Vb-ex

1 3

Controlling transient gas flow inreal‑world pipeline…

nodes

Vw1

has to be smaller than or equal to the flow over the second set of boundary

nodes

Vw2

. These conditions must be met if an associated flow direction

f=(f+

,f−)

is active. By definition the node sets

Vw1

and

Vw2

have to be a subset of either

f−

, while they do not need to belong to the same set.

An example of a relation requiring the constraints is as follows: Consider a trans-

port pipeline intersecting a network station. These constraints can then be used to

model a flow direction that not only restricts the flow at the boundary nodes, but

also enforces that gas is directed into the intersecting transport pipeline. We do this

by stating that the outflow of the transport pipeline at the network station’s boundary

is larger than its inflow into the network station.

Note that the flow over a set of boundary nodes

for time t is defined as

∑v∈V

�

v,t�

, which is potentially a nonlinear expression due to the absolute value.

However, by the above definition, the boundary node sets

Vw1

and

Vw2

of some

(

−

),V

1,V

)

∈

are subsets of either

f−

. For this reason, we always

know the sign of the flow over

in advance, and hence using the following func-

tion definition for easier notation

we can write the constraints for each

t∈T

as:

2.8 Operation mode transition times

In addition to the characterization of operation modes given previously, there exists

additional restrictions, namely those concerning transition times between any two

operation modes. In contrast to the rest of Sect.2, we only give a verbal descrip-

tion and do not present a MIP model formulation here, since none of the MIP vari-

ants solved in our heuristic algorithm uses these constraints, see Sect.4. Instead, it

explicitly checks potential solutions for complying with this restriction.

Changing the operation mode of a network station can be a complex procedure

if the two involved operation modes differ significantly in terms of the actively used

elements. The changes in terms of modes and configurations that are needed for the

transition are often required to happen in a particular order. Furthermore, it may be

necessary to switch into a third, intermediate operation mode during the change to

ensure safe operation at all times. However, for the scope of this paper, we are only

interested in the implications of the transition process. Therefore, we define a feasi-

ble transition from operation mode A to mode B as follows: We allow the transition

to happen instantly at some time point

t∈T0

. However, we consider the time

𝜃(A,B)

of the transition process from A to B by demanding that the network station is in

operation mode A in the time interval

[

t−

𝜃(A,B)

)

preceding t and in mode B for

sgn

∶F×Vb→{−1, 1},

(

(f+,f−),v

)

→

{

1 if v∈f

−1 if v∈f−

(38)

f,t=1→

∑

v∈Vw1

sgn(f,v)dv,t≤

∑

v∈Vw2

sgn(f,v)dv,t

∀(f,V

)∈W.

F.Hennings et al.

1 3

the subsequent time interval

[

t,t+

𝜃(A,B)

. In addition, this whole transition period

[

t−

𝜃(A,B)

,t+

𝜃(A,B)

is reserved for this transition to happen, i.e., two transition peri-

ods cannot overlap. Note that we allow the last transition period to exceed our time

horizon, which makes operation mode changes up to the very end of the time hori-

zon possible.

As an example, we denote in Fig. 2 an operation mode sequence over time

that includes the involved transition times. Here, each time point

represents the

switching time of the corresponding transition into the operation mode X, i.e., the

first

t∈T0

in which operation mode X is active. Note that the displayed example

includes a conflict in the transitions from mode C to mode D and from mode D to

mode E since the two corresponding transition periods overlap. Hence this sequence

of operation modes does not represent a feasible solution.

2.9 Scenario andinitial state

For the future, we are given scenario values for the station’s boundary nodes in

terms of pressure and inflow (where outflow is represented as negative inflow, see

the beginning of Sect.2). While we have one pressure value

pv,t

per boundary node

v∈Vb

for each future time point

t∈T

, the flow demands are only given for sets of

boundary nodes. These sets are called the fence groups of the network station and

form the set

. For each

g∈G

and time point

t∈T

, the sum of the inflow values of

all boundary nodes in g should be equal to the demand value



dg,t

We do not require strict obedience of the given demand values

pv,t

and



dg,t

, but

instead allow deviations from them, which we punish in the objective function.

These deviations are captured in the slack variables

𝜎

. For pressure, we have

𝜎p+

v,t

and

𝜎p−

v,t

capturing the positive and negative difference between the pressure values

of boundary node v at time t and the given demand

pv,t

. The slack variables

𝜎d+

v,t

and

𝜎d−

v,t

are associated with the inflow. They capture the positive and negative contribu-

tion to the difference between the inflow demand



dg,t

of fence group

g∈G

at time

step t and the sum of the corresponding inflow variables of each boundary node v of

the fence group g.

We model the described relations for each future time step

t∈T

as:

Time

tBtCtDtE

AB CDE

A→B

B→C

C→D

D→E

Fig. 2 Transition time example with five operation modes and four transitions, from which two conflict.

The time point

represents the switching time of the corresponding transition into the operation mode

X, i.e., the first

t∈T0

in which operation mode X is active

1 3

Controlling transient gas flow inreal‑world pipeline…

Note that by using inflow slacks, it is possible to choose a flow direction that does

not fit the inflow demand values



dg,t

in terms of the Eqs.(34–35).

The second set of prescribed values are those of the initial state represented by time

0. All variables for this time are actually parameters of the model and fixed to the cor-

responding value.

2.10 Objective

The objective function covers two different aspects. On the one hand, it should mini-

mize the deviation from the given future scenario in terms of inflow and pressure at

the boundary nodes. On the other hand, we favor those solutions with stable control

of the single elements and hence try to reduce the number of needed control changes.

The latter is preferable for a couple of reasons. First, a limited number of changes ena-

bles the gas network operators to understand, evaluate, and actually perform the desired

control recommendations. Due to the inertia of the gas, they are used to operating the

network using a small set of essential adjustments rather than dynamically performing

small changes all the time. Finally, a stable control also reduces strain on the technical

elements.

To define a measure for the stability, we first quantify the discrete changes of the

control in binary variables. The variable

𝛿om

captures the change of the network station

into a new operation mode at time t, while the variable

𝛿rg

a,t

represents the change into a

new mode of regulator a at time t. We have to track the regulator mode changes explic-

itly, since their modes are not prescribed by the operation modes, see the description of

the network station structure in Sect.2.7 for an explanation. Furthermore, we capture

the start-up of compressor unit u at time t in the variable

𝛿us

u,t

. Since starting a compres-

sor is a very time and energy intensive action, it should be avoided if possible. The fol-

lowing set of constraints guarantees for each future time step

t∈T

that the variables

attain the value 1 if a corresponding change occurs and 0 otherwise.

(39)

p

v,t

−𝜎

v,t

+𝜎

p−

v,t

∀v∈V

(40)



g,t=

∑

v∈g(

dv,t−𝜎d+

v,t+𝜎d−

v,t

)

∀g∈G

(41)

𝛿om

≥

mom

o,t−mom

o,t−1∀o∈O

(42)

𝛿om

≤2

−

o,t−

o,t−1∀

∈O

(43)

𝛿rg

a,t≥m

a,t

−m

a,t−1

∀a∈A

(44)

𝛿rg

a,t

≤2−m

a,t

−m

a,t−1

∀a∈A

F.Hennings et al.

1 3

In Eq.(49), the set

Cu,a

denotes the subset of configurations

of the compressor sta-

tion

a∈Acs

, for which compressor unit

u∈Ua

is turned on.

In order to obtain a smooth operation for network situations without discrete

mode switching, we add variables tracking the change of the operation point of sin-

gle elements, i.e., their corresponding changes in flow, incoming pressure and out-

going pressure. We do this for all elements that can precisely control the point of

operation, i.e., regulators and compressor stations in active mode, and for all times

without discrete mode switching, i.e. a stable network station operation mode or reg-

ulator mode. The variables

𝛿rg-pl

a,t

𝛿rg-pr

a,t

, and

𝛿rg-q

a,t

represent changes of the incoming

pressure, outgoing pressure, and flow of an active regulator a, while

𝛿cs-pl

a,t

𝛿cs-pr

a,t

, and

𝛿cs-q

a,t

represent the corresponding value changes of an active compressor station a.

We establish the behavior of these variables using the following constraints for each

(

𝓁

,r)=a∈Arg ∪Acs

and each

t∈T

(45)

𝛿rg

a,t≥m

a,t−m

a,t−1

∀a∈A

(46)

𝛿rg

a,t

≤2−m

a,t

−m

a,t−1

∀a∈A

(47)

𝛿rg

a,t

≥m

a,t

−m

a,t−1

∀a∈A

(48)

𝛿rg

a,t

≤2−m

a,t

−m

a,t−1

∀a∈A

(49)

𝛿us

u,t≥

∑

c∈Cu,a

mcf

c,a,t−

∑

c∈Cu,a

mcf

c,a,t−1∀u∈Ua∀a∈A

𝛿

t∈{0, 1}

𝛿rg

a,t∈{0, 1}∀a∈Arg

𝛿us

u,t

∈{0, 1}∀u∈U

∀a∈Acs

(50)

𝓁,t−p𝓁,t−1≤𝛿

X-pl

a,t+(m

a,t+mcl

a,t+𝛿)(p𝓁,t−p𝓁

,t−1)

(51)

𝓁,t−1−p𝓁,t≤𝛿

X-pl

a,t+(m

a,t+mcl

a,t+𝛿)(p𝓁,t−1−p𝓁

)

(52)

pr,t−pr,t−1≤𝛿

X-pr

a,t+(m

a,t+mcl

a,t+𝛿)(pr,t−p

r,t−1)

(53)

pr,t−1−pr,t≤𝛿

X-pr

a,t+(m

a,t+mcl

a,t+𝛿)(pr,t−1−p

r,t

)

1 3

Controlling transient gas flow inreal‑world pipeline…

In these constraints, we define X=rg as well as

𝛿

=𝛿

a,t

for

a∈Arg

, and X=cs as

well as

𝛿=𝛿om

for

a∈Acs

Note that we needed to add the upper bound constraints (42), (44), (46), and (48)

for the discrete change variables

𝛿om

and

𝛿rg

a,t

to allow them to be equal to 1 if and

only if there is a discrete change. Otherwise, it would be possible to set them to

1 independent of the corresponding mode variable values, thereby deactivating the

associated constraints (50)–(55) and avoiding potentially higher costs imposed by

the continuous change variables.

Finally, we can state our objective function, which minimizes the weighted sum

of the change variables as well as the slack variables defined in Sect.2.9 as

Here, the

w∗

parameters denote the corresponding positive weights given to the sin-

gle quantities. Note that the weights for pressure and flow slack variables are addi-

tionally multiplied by the length of the corresponding time interval to represent the

amount of deviation over time correctly.

2.11 Final model

Putting everything together, we can formulate our problem in the following transient

gas flow model

(54)

a,t−qa,t−1≤𝛿

X-q

a,t+(m

a,t+mcl

a,t+𝛿)(qa,t−q

a,t−1)

(55)

a,t−1−qa,t≤𝛿

X-q

a,t+(m

a,t+mcl

a,t+𝛿)(qa,t−1−q

a,t

)

(56)

min ∑

t∈T

((

𝜏(t)−𝜏(t−1)

)∑

v∈Vb

(

w𝜎-p ⋅(𝜎p+

v,t+𝜎p−

v,t)+w𝜎-d ⋅(𝜎d+

v,t+𝜎d−

v,t)

)

+wom ⋅𝛿om

t+∑

a∈Arg

wrg ⋅𝛿rg

a,t+∑

u∈Ua,a∈Acs

wus ⋅𝛿us

u,t

+∑

a∈Arg (wrg-pl ⋅𝛿rg-pl

a,t+wrg-pr ⋅𝛿rg-pr

a,t+wrg-q ⋅𝛿rg-q

a,t)

∑

a∈A

(

wcs-pl ⋅𝛿cs-pl

a,t+wcs-pr ⋅𝛿cs-pr

a,t+wcs-q ⋅𝛿cs-q

a,t

))

F.Hennings et al.

1 3

Note that we apply the constraints only starting from time step 1, explicitly exclud-

ing the initial time step 0. We do this since the initial pressure values, flow values,

and modes described in Sect.2.9 are not guaranteed to fit our model, and only serve

as a starting point for the calculations.

3 Feasible operating range ofcompressor station configurations

In the following section, we describe the creation of the feasible operating range

polytope for each configuration

c∈Ca

of a compressor station

a=(𝓁,r)

. This poly-

tope is represented as the intersection of half-spaces

and was already used in the

model description for compressor stations in Sect.2.1 in Eq.(12).

min (56)

s.t.

∀t∈T(1)−(4),(8)−(11),(30)−(31),(50)−(55)∀a∈A

(5)−(7),(32)∀a∈Acs ∀c∈Ca

(12)∀a∈Acs ∀c∈Ca∀(𝛼0,𝛼1,𝛼2,𝛼3)∈Hc

(13)−(15)∀a∈Api

(16)∀a∈Ars

(17)−(20),(29)∀a∈Ava

(21)−(25),(43)−(48),(50)−(55)∀a∈Arg

(26),(34)−(35),(39)∀v∈Vb

(27)∀v∈V0

(28),(33)

(36),(41)−(42)∀o∈O

(37)∀v∈Vb-ex

(38)∀w∈W

(40)∀g∈G

(49)∀a∈Acs u∈Ua

{0, 1}∋mby

a,t,mcl

a,t∀a∈Acs,

mcf

c,a,t∀a∈Acs ∀c∈Ca,

mop

a,t∀a∈Ava,

mcl

a,t,mby

a,t,mac

a,t,𝛿rg

a,t∀a∈Arg,

mom

o,t∀o∈O,

mfd

f,t∀f∈F,

𝛿om

𝛿us

u,t

∀u∈U

∀a∈Acs.

1 3

Controlling transient gas flow inreal‑world pipeline…

We shortly recall the structure of a compressor station arc from Sect.2.1. Each

compressor station

a∈Acs

is associated with a set of compressor units

, which

each represent the union of a compressor machine and the power providing drive.

These compressor units can be combined in series and/or parallel to either allow

for a higher compression ratio, a higher flow rate, or a mixture of both. The set of

all available serial-parallel combinations is called the set of configurations

of a

compressor station.

To create the set

for

c∈Ca

, we first obtain a feasible operating range polytope

for each of the compressor station’s compressor units. This polytope includes a lin-

ear approximation of the nonlinear constraint resulting from the maximum power

bound of the unit’s drive. Afterwards, we sketch the procedure of Walther and Hiller

(2017) to combine these polytopes and obtain

Note that we drop the index a for the compressor station as well as the time index

in this section for ease of notation.

3.1 Feasible operating range forasingle compressor unit

For each compressor unit, we are given a feasible operating range of the correspond-

ing compressor machine as a polytope in the space

(

p𝓁

)

, where the volumetric

flow rate

is defined as

Note that we consider the density

𝜌

and the pressure

at the incoming node

𝓁

of the

compressor. For this reason, we use for compressor units the incoming node’s com-

pressibility factor

𝓁

∶= z(p

𝓁

,0)

instead of the average value

used for pipes and

resistors as defined in Sect.2.2.

Usually, the feasible operating range, sometimes also referred to as “character-

istic diagram” or “performance curve”, is given as area in the dimensions

(Q,Had)

restricted by a set of possibly concave quadratic curves, see ,e.g., the description of

Fügenschuh etal. (2015) or Odom and Muster (2009). The quantity

Had

denotes the

specific change in adiabatic enthalpy and is defined as

using for the isentropic exponent

𝜅

the constant value 1.296, as stated by Fügens-

chuh etal. (2015). We can transform such an operating range into our format by

applying the unique transformation from

Had

𝓁

obtainable from (57).

In addition to the feasible operating range polytope, each compressor unit is given

an upper bound on the absolute pressure increase



𝛥p

≥p

−p

𝓁

achievable by the

compressor machine and an upper bound on the maximum power



that may be

delivered by the drive. The power needed for compression is defined as

=q∕𝜌𝓁, with 𝜌𝓁=

𝓁

RsTz𝓁

(57)

ad =RsTz𝓁

𝜅

𝜅−1

[(

p𝓁

)

𝜅−1

𝜅

−1

F.Hennings et al.

1 3

Here,

𝜂ad

denotes the adiabatic efficiency of the compression, which we assume to

be a known constant per compressor unit. We give an example of a feasible operat-

ing range of a compressor unit combined with the two constraints in Fig.3a, where

different levels of the maximum pressure difference bound and the maximum power

constraint are denoted based on different values for the incoming pressure

p𝓁

We now create a first version of the feasible operating range of the compressor unit

in the space

(p𝓁,pr,q)

that does not include the maximum power constraint. Therefore,

we transform each of the faces

𝛼

0+𝛼1Q+𝛼2

p𝓁

≤

of the original polytope in the

space

(

𝓁

)

into constraints

𝛼 0p𝓁+𝛼 1pr+𝛼 2q

≤

of the desired higher-dimensional

space:

To bound the polyhedron described by the new constraints, we add the restriction of

the absolute pressure difference as well as two pressure bounds of the end nodes of

the compressor station arc

(𝓁,r)=a∈Acs

containing the compressor unit.

(58)

qHad

𝜂ad

RsTz𝓁

𝜅

𝜅−1

[(

p𝓁

)𝜅−1

𝜅

−1

𝛼0+𝛼1Q+𝛼2

p𝓁

≤

⇔

𝛼0+𝛼1

qRsTz𝓁

p𝓁

+𝛼2

p𝓁

≤

⇔

𝛼0p𝓁+𝛼1RsTz𝓁q+𝛼2pr≤

⇔

𝛼

0p𝓁+

𝛼

1pr+

𝛼

≤

Fig. 3 The feasible operating range of a compressor unit. The grey region shows the operating range

given as a polytope. For different values of the incoming pressure

p𝓁

the blue dashed lines represent the

upper bound on the absolute pressure increase



𝛥p

, and the red solid lines illustrate the power bound



While the left picture shows the original nonlinear nonconvex constraint resulting from the power bound,

the right picture shows the linearized version, see Sect.3.2

1 3

Controlling transient gas flow inreal‑world pipeline…

A picture of the corresponding three-dimensional polytope created based on the fea-

sible operating range of Fig.3a is given in Fig.4a. Note that it does not include the

max power constraint yet.

3.2 Power bound linearization

To complete the three-dimensional feasible operation range polytope, we now add

the maximum power constraint described in 3.1. Note that this constraint

P≤



cuts

into the original two-dimensional feasible operating range in a nonlinear and non-

convex fashion, see Fig.3a. Since the same holds for the feasible operating range

representation in

(p𝓁,pr,q)

, we derive a linear approximation to this constraint.

For the approximation, we generate a set of N random sample points from within

the three-dimensional operating range polytope of the compressor unit derived in

the previous section. For each point, we determine the corresponding compression

power using Eq.(58). Given this, we apply an ordinary least-squares method in

order to determine the coefficients of the linear representation of

. This enables us

to ultimately add the constraint in the form

to the three-dimensional polytope of the operating range and thereby complete its

definition. An example of the final polytope is illustrated in Fig.4b, while the lin-

earized power bound projected to the original two-dimensional operating range can

be seen in Fig.3b.

−p𝓁≤



𝛥pp𝓁≥p

𝓁

pr≤p

𝛽0

+𝛽

𝓁

+𝛽

q≈P≤



P𝛽

,𝛽

∈

ℝ

Fig. 4 The feasible operating range of a compressor unit in the space

(p𝓁,pr,q)

, computed from the two-

dimensional operating range shown in Fig.3. While the left picture shows the lifted polytope based on

the original feasible operating range, the maximum absolute pressure difference, and the end nodes pres-

sure bounds, the right picture also includes the linearized constraint resulting from the power bound, see

Sect.3.2

F.Hennings et al.

1 3

3.3 Feasible operating range foracompressor station configuration

As final step, we derive the polytope description for each configuration by combin-

ing the polytopes of the used compressor units. Each configuration is given as a

serial sequence of parallel compressor machine arrangements. Therefore, we first

combine the compressor unit polytopes of each parallel arrangement and then cre-

ate the final configuration polytope from these. The set of facets of this configura-

tion polytope ultimately defines the set of half-spaces

. The overall procedure was

introduced by Humpola etal. (2015), while we exactly followed the steps of the

variant presented by Walther and Hiller (2017).

4 Specialized network station algorithm

The MIP model

presented in Sect.2.11 turns out to be quite challenging accord-

ing to our experiments, despite only considering a part of a larger gas network by

restricting ourselves to a single network station. For this reason, we create a special-

ized heuristic algorithm to solve the problem. It should be able to efficiently solve

the problem even on harder instances, i.e., those featuring complex network stations

and many discrete time steps. Furthermore, we require it to reliably find feasible

solutions for real-world instances and achieve a solution quality that is close to the

optimum.

To reach these goals, we design the algorithm based on the before mentioned

observation that pipes inside a network station are relatively short due to the proxim-

ity of the contained elements. This affects the pipes’ capability to store gas, which is

often referred to as linepack. The consideration of linepack is the main improvement

of a transient pipe model over a stationary one. However, for short pipes, linepack is

insignificant in comparison to the long pipelines outside of network stations, making

the presented transient optimization problem a “stationary-flavored” one.

This property inspired us to determine the network station’s operation modes and

flow directions based on the solution of stationary (or steady-state) variants of

solved

independently for each time step. The above description implies that these values

should be similar to the ones of the optimal solution of the overall transient problem.

Furthermore, the temporal decoupling results in a significant complexity reduction of

the corresponding MIP models. In the algorithm, we ensure that the sequence of opera-

tion modes over time respects the transition time conditions introduced in Sect.2.8,

which we explicitly did not include in

After determining the operation modes and flow directions, we obtain our final solu-

tion by using another variant of

, which covers multiple time steps. By prescribing

the values determined in the stationary calculations and solving the model in a rolling

horizon fashion, we make this variant tractable as well. This procedure is called tran-

sientSmoothing and is the topic of Sect.4.3. The whole algorithm is summarized

in Algorithm1, in which we split the creation of the sequence of operation modes and

flow directions into the two-step procedure that we will present in Sect.4.2.

1 3

Controlling transient gas flow inreal‑world pipeline…

Note that our algorithm has no guarantee of global optimality. We step-wise deter-

mine the decision variables without considering all potential consequences for the

subsequent routines, which may lead to sub-optimal solutions. However, the order in

which these decisions are made matters. Algorithm1 first determines the operation

modes and flow directions in the subroutines initialSolutionCreation and

improvementHeuristic, while taking into account the deviation from the given

future demands. The objective function weights associated with these decisions are

wom

for changing operation modes,

wus

for starting new compressor units, as well as

w𝜎-p

and

w𝜎-d

for not fulfilling the demand requirements. Usually, these weights are the most

dominant ones, see the used objective function values in our computational experi-

ments stated in Sect.5.1. The remaining parts of the solution are obtained in the subse-

quent step transientSmoothing. Since here the operation modes and flow direc-

tions are already fixed, the remaining decisions may lead to higher values in the penalty

variables, like for example those tracking the continuous change in the points of opera-

tion. However, due to the smaller objective weights these variables also have a smaller

influence on the overall objective value. Hence we can summarize that Algorithm1 is

likely to obtain good overall solutions for the given problem, which is confirmed by the

computational results we present in Sect.5.

4.1 Model variants

As mentioned above, we do not solve the transient model

directly, but rather use

the following variants in our overall solution approach.

4.1.1

—Stationary model

For the first variant, we solve a stationary version of the model to determine the

best operation mode and flow direction for one independent time step t. The neces-

sary changes in

mainly affect the pipes. Here, Eq.(13) is no longer part of the

model. Furthermore, a pipe can no longer store gas, which balances the incoming

and outgoing pipe flows, i.e.,

q𝓁,a,t=qr,a,t=qa,t

holds for all pipes

(

𝓁

,r)=a∈

Thus, we can state the stationary version of the Momentum Eq. (15) for pipe

(

𝓁

,r)=a∈

and stationary time step t as

For all other elements, the constraints only consider variables of the same time step.

Therefore we simply apply these for the time step t of the stationary case. As a final

step, we adjust the objective function. We keep the slack penalties for each boundary

node, which are defined based on individual time steps. In addition, we still track the

r,t−p𝓁,t+

𝜆

v𝓁,a

vr,a

qa,t+

(

p𝓁,t+pr,t

)

F.Hennings et al.

1 3

change to a new network station operation mode in variable

𝛿om

as well as the start

of new compressor units using

𝛿us

u,t

for all

∈U

a∈A

. This is done by calling the

model with the parameter prevMode, which represents the operation mode of the

previous time step, see Algorithm2. All the other variables that track changes are

not considered. Therefore, we have the following stationary objective function for

the time step t:

4.1.2

—Transient withfixed operation modes andflow directions

In this variant, we are given for all future time steps

t∈T

a fixed operation mode

and flow direction

of the network station. This removes the corresponding varia-

bles

mom

ot,t

and

mfd

ft,t

from the model by setting them to 1 or 0, depending on whether

they are chosen or not. Since the modes of valves and compressor stations, as well as

the configurations of active compressor stations, depend on the selected operation

mode of the network station, their corresponding variables are eliminated as well.

Furthermore, all indicator constraints that use the binary variables mentioned above

in their logical condition are resolved, i.e., we either remove them from the model if

the condition is not fulfilled or add the implied constraint otherwise. Finally, we no

longer need to use a disjunctive model for compressor stations. Instead, we apply the

corresponding constraints(2)–(12) for the chosen mode and configuration of time

step t directly to the variables

𝓁

pr,t

, and

qa,t

of each compressor station

(

𝓁

,r)=a∈Acs

4.1.3

Psf

—Stationary withfixed operation mode

For this last version of model

, we combine the two variants above and use the

stationary version of the model for one independent time step t with an already fixed

operation mode. Note that in contrast to model

, the flow direction of the network

station is not determined yet, which results in more binary variables and more indi-

cator constraints that are not yet resolved.

min (

𝜏(t)−𝜏(t−1)

)

⋅

∑

v∈Vb

w𝜎-p ⋅(𝜎

v,t+𝜎

p−

v,t)+w𝜎-d ⋅(𝜎d+

v,t+𝜎d−

v,t

)

+wom ⋅𝛿om

∑

u∈U

,a∈Acs

wus ⋅𝛿us

u,t.

1 3

Controlling transient gas flow inreal‑world pipeline…

4.2 Determining operation modes andflow directions

In this section, we present the procedure that finds a series of operation modes of

the network station over time by solving the stationary variants of

. We split the

approach into two steps. First, we create an initial solution by a greedy, forward-

oriented procedure presented in Sect.4.2.1. It represents the function called ini-

tialSolutionCreation in Algorithm1. In a second step, we try to improve

this sequence of operation modes. Here, we replace certain operation modes of

the sequence by similar ones in order to obtain a better solution. This second step

is referred to as improvementHeuristic in Algorithm1 and is described in

Sect.4.2.2. For the sake of efficiency, we try to reduce the number of times we need

to solve model

and try to use model

Psf

instead. In

Psf

, we prescribe an opera-

tion mode for the network station, which makes the remaining problem very easy

and extremely fast to solve. In contrast to this, the stationary model

considering

all available operation modes is solved significantly slower, especially when used on

larger network stations. Therefore, we only use it if necessary.

Besides the operation mode, we return for each time step t a flow direction. The

returned flow direction for t is the one determined by the stationary model that com-

puted the optimal operation mode for t. For simplicity of notation, we omit the flow

directions in the returned entities in the remainder of this section.

4.2.1 Initial solution creation

To find a first feasible sequence of operation modes over time, we follow a rather

simple idea: We determine an operation mode for time step t by first testing the used

operation mode of the previous time step

t−1

using

Psf

. The general stationary

model

to determine the best operation mode for t is only called when the previ-

ously used operation mode is infeasible or yields a costly solution in terms of the

objective function value. Using this approach, we keep the number of needed opera-

tion mode changes small and speed up the algorithm by reducing the number of

calls to

. A detailed description is given as Algorithm2.

F.Hennings et al.

1 3

In Algorithm2, we call

with the parameter validModes, which replaces the

set of valid network station operation modes

. Furthermore, we call the functions

modeAvailable, transitionsWork, and notSoonInfeasible in Algo-

rithm2. The function modeAvailable deals with the operation mode unavail-

ability introduced in Sect.2.7, and checks if a given operation mode o is available

at time t. In transitionsWork a given operation mode sequence is tested for

validity regarding the transition time restriction presented in Sect.2.8. Finally, we

check in notSoonInfeasible, if choosing a new operation mode o for time t

would result in an infeasibility at one of the subsequent time steps that is caused by

a combination of operation mode unavailability and too long transition times. More

specifically, we check if the new operation mode o for time t will become unavail-

able in one of the future time steps. If this is the case, we test if there is enough time

left to transition into another operation mode until then. This also take into account

the time needed by the transition from the old operation mode at time

t−1

to this

new operation mode o at time t.

Note that Algorithm2 may abort without a feasible solution at line14. However,

this is extremely unlikely to happen in real-world instances. A possible cause for

failing at this line is that the algorithm may, due to a restrictive combination of oper-

ation mode availability and transition times, not find a valid operation mode to add

1 3

Controlling transient gas flow inreal‑world pipeline…

to the set validModes before the call of the model

. There are multiple rea-

sons why we do not expect this situation to happen. Firstly, the compressor units are

rarely unavailable in practice, and even when all units of a compressor station are

unavailable, the compressor station is still operational in bypass or closed modes.

Secondly, the restrictions imposed by the transition times are usually rather mild in

the sense that long transition times only occur for select pairs of operation modes.

Finally, the function notSoonInfeasible prevents choosing an operation mode

that may lead to a conflict in upcoming time steps. Hence it is likely that valid-

Modes contains at least one operation mode. A further possible cause for abort-

ing the algorithm occurs when

is infeasible for all given operation modes. This

again is not expected for real-world instances, namely because the demand scenario

in terms of pressure and flow is adjustable without limitation, even if these adjust-

ments would result in a massive usage of the corresponding slack variables. Hence,

the model can only be infeasible if all available operation modes have a conflict in

every possible demand scenario. Therefore, infeasibility is a sign of a structural con-

flict in the network elements’ modes and configurations imposed by the correspond-

ing operation modes of the network station and is not assumed to happen for well-

formed, realistic operation modes. The fact that Algorithm2 has never aborted in

one of our real-world test cases presented in Sect.5.3 supports our argumentation.

However, since we cannot guarantee feasibility here, we leave the removal of this

theoretical weakness of Algorithm2 open for future research.

F.Hennings et al.

1 3

4.2.2 An improvement heuristic

After having found an initial feasible solution, we look for further improvements

using Algorithm3. Therefore, we identify all sequences of identical operation

modes over time in the solution. We call these sequences stable phases, or just

phases, of a feasible solution. Obviously, the switch from one phase to the subse-

quent one happens when changing the operation mode. For each of these phases,

we check if its operation mode can be replaced by a similar one that is more ben-

eficial in terms of the objective function value. The general idea is to compare

different network station operation modes over longer periods of time, which is in

contrast to the focus on individual time steps implemented in Algorithm2. Note

that we only use the fast solving MIP model variant

Psf

here.

To generate new operation modes to check for improvements, we call the

function convexCombination, which is the main feature of Algorithm 3.

To define it, we use the function M(o,a) that returns the mode or active con-

figuration of a valve or compressor station a prescribed by operation mode o,

see Sect.2.7. Furthermore, we denote by

U(x)

the compressor units used in mode

or configuration

x∈{by, cl}∪Ca

for some compressor station

a∈Acs

, where

U(by)=U(cl)=�

. We first define the function convexCombinationCS

on a

tuple (x,y) with

x,y∈{by, cl}∪Ca

Note that convexCombinationCS

a(x,y)⊆{by, cl}∪Ca

holds. Now we can

define convexCombination on a tuple

(o1,o2)

of operation modes as

For a compressor station a, we permit it to be in a configuration using a compres-

sor unit set “in-between” the used compressor unit set of the configurations used in

and

for a. However, we only allow the exact valve mode combination used in

or the one used in

. This choice is based on the structure of operation modes

described in Sect.2.7. While changing the configuration of an active compressor

station in an existing operation mode usually yields another existing operation mode,

the set of valid combinations of valve modes is rather limited. Therefore, finding a

valid valve setting from the combination of valve modes of two operation modes is

unlikely and hence not included in the algorithm.

Algorithm3 uses the function transitionsWork, which works in the same

way as described above for Algorithm2. Furthermore, it calls allModesAvail-

able, which is similar to modeAvailable from Algorithm2. Here, we check the

availability of a whole sequence of operation modes at those times corresponding

𝚌𝚘𝚗𝚟𝚎𝚡𝙲𝚘𝚖𝚋𝚒𝚗𝚊𝚝𝚒𝚘𝚗𝙲𝚂

(x,y) ∶= {x,y}∪

{

c∈Ca

∀u∈U(x)∩ U(y)∶u∈U(c)

)

∧

(

∀u∈U(c)∶u∈U(x)∪U(y)

)}.

𝚌𝚘𝚗𝚟𝚎𝚡𝙲𝚘𝚖𝚋𝚒𝚗𝚊𝚝𝚒𝚘𝚗

(o1,o2) ∶=

{

o∈O|

(∀a∈Ava ∶M(o,a)=M(o1,a)∨∀a∈Ava ∶M(o,a)=M(o2,a))

∧∀a∈Acs ∶M(o,a)∈𝚌𝚘𝚗𝚟𝚎𝚡𝙲𝚘𝚖𝚋𝚒𝚗𝚊𝚝𝚒𝚘𝚗𝙲𝚂a

(

M(o1,a),M(o2,a)

)}.

1 3

Controlling transient gas flow inreal‑world pipeline…

to their positions in the sequence. Furthermore, we evaluate the objective function

value of a sequence of operation modes using the function obj by successively

calling

Psf

for each operation mode and time corresponding to its position in the

sequence. If one of the models turns out to be infeasible, the returned objective

function value is infinity.

Finally, we start the algorithm in the backward oriented mode since the initial

solution is obtained by Algorithm2, which works in the forward direction. Further-

more, we highlight that Algorithm3 has the potential to reduce the total number of

needed operation mode changes since the two original operation modes defining the

change from one phase to the other are always part of the result of convexCombi-

nation. Since it is not possible to introduce new operation mode changes, we can

even guarantee the new number of operation mode changes to be less than or equal

to that in the initial solution.

4.3 Transient solution smoothing

As a final step of our specialized network station approach presented in Algorithm1,

we solve the transient model variant

with fixed network station operation modes

and flow directions. We obtain the operations modes as well as the flow directions

for each time step from the stationary model solutions created in the previous steps

of the algorithm.

In our computational experiments, we observed that, even though most of the

binary decision variables of

are fixed in

, only a limited number of time

steps can be solved for large network stations. Therefore, we use a rolling horizon

approach to solve

, which we describe in Algorithm4. The idea is to solve a

series of smaller

models with a fixed time horizon size h representing the

F.Hennings et al.

1 3

number of time steps for which we solve, which include the time step for the given

initial state. After each iteration, we save the solution regarding the earliest time step

for the final result and shift the time horizon by 1 to prepare for the next call of

In the function call to solve

in Algorithm4, we pass the operation modes and

flow directions that correspond to the current time horizon as well as the fixed initial

state.

The main benefit of this method is that increasing the size

|T0∶= {0, …,k}|

of the overall time horizon only increases the number of equally sized MIP mod-

els to solve, which hence have similar complexity. Therefore, the increase of

|T0∶= {0, …,k}|

is expected to result in a proportional run time increase. If we

would have solved

over the entire time horizon, increasing the time horizon

would result in a larger complexity of the model, which may lead to an exponential

increase in runtime instead.

5 Computational experiments

When designing Algorithm1, we wanted it to be competitive in terms of the reli-

ability to find a solution, the quality of that solution, and the overall execution time.

In order to verify that we indeed achieved these objectives, we evaluated the algo-

rithm on a large number of test instances. These instances represent network stations

in the network of our project partner Open Grid Europe(OGE), for which we gener-

ated scenario values based on historic real-world situations.

5.1 Instances

We considered 7 different network stations from the network of OGE with different

sizes and properties. An overview can be found in Table1.

For each of these stations, we have a set of 159 instances. To create these

instances, we solved a macroscopic gas flow model on an aggregated version of

Table 1 Overview of different properties of the 7 network stations A to G

For station G, the topology changed during the considered time period. Hence, we denoted by G

min

the

minimal values of all the quantities for station G and by G

max

the corresponding maximum values

Name

|V|

|A|

∑

∈Api

�

Api

�

(km)

|Ca|

∀a∈

Acs

|O|

|F|

b-ex|

|W|

A 14 11 0.001 16 34 3 1 0

B 11 12 0.001 18 23 2 1 0

C 27 34 0.012 2 13 4 2 0

D 25 31 0.404 2, 6 92 6 2 0

E 48 67 0.308 3, 5 82 12 1 4

F 51 66 0.024 2, 3, 6, 12 2836 3 3 0

min

118 148 0.079 1, 1, 2, 7, 20 1267 15 2 0

max

120 150 0.079 1, 1, 2, 7, 20 1285 20 2 0

1 3

Controlling transient gas flow inreal‑world pipeline…

the complete network using simplified versions for all the network stations. A full

description of this model is given by Hoppmann etal. (2019).

Each of the 159 instances corresponds to one point in time, for which we are

given the actually occurred state of the network of our project partner. These points

in time are distributed over 3 months and have a time difference of at least 12 hours.

In addition to one of these states, we have measured pressure and flow values for the

following 12 hours at the boundary nodes of the entire network. The macroscopic

gas flow model can then be formulated using the state of the network as initial state

and the pressure and flow values for the following 12 hours as demand values at the

boundaries of the network. The solutions of the macroscopic gas flow model for

each of these instances yield pressure and inflow values at the boundaries of each

network station. The input data of one instance in the test set of this paper hence

consists of a) the topology of one of the 7 network stations, b) the initial state values

for the elements contained in this station for one of the 159 given points in time, and

c) the pressure and inflow values at the network station’s boundaries for the follow-

ing 12 hours taken from the corresponding macroscopic gas flow model solution.

For the macroscopic gas flow model, the time horizon of 12 hours consists of 4

periods of 15 minutes followed by 11 periods of 60 minutes, resulting in 15 future

time steps in total. For our analysis, we created the following four different partitions

of the 12 hour time horizon built from left to right:

Our partitioning of the time horizon can have additional or missing time points com-

pared to the 15 time steps used in the macroscopic model. In order to create scenar-

ios values that fit the time horizon, we interpolated using the original values.

To conclude, we used the following set of weights for the objective function for

all instances for X

∈{rg, cs}

The values regarding flow are based on the unit

1000m3∕h

for volumetric flow under

normal conditions. Using the linear transformation

[

1000m3∕h]=

3600

1000𝜌

0[kg∕s

]

based on the given normal density

𝜌0

parameter, we can transform the corresponding

values into mass flow. The unit of

w𝜎-d

is deduced by

1000 m

∕h

⋅

1000 m

To give an impression of the instances’ complexity, we list in Table2 the size of

the original model

for each network station and each of the above-defined time

horizon partitions. The given numbers state the number of continuous variables,

binary variables, and constraints, each rounded to the nearest hundred. As one can

see, the values increase proportionally to the number of time steps. Note that we did

not solve

directly but rather used our Algorithm1 presented in Sect.4 to obtain

12 time steps: 4 ×15 min, 5 ×60 min, 3 ×120 min

24 time steps: 4

×15 min, 18 ×30 min, 2 ×

60 min

48 time steps: 48

×15 min

96 time steps: 96 ×7.5 min.

w𝜎-p

=1000.0∕(bar h)w

𝜎-d

=100.0∕(1000 m

=1000.0 w

1200.0

wrg

=50.0

wX-pl

wX-pr

=10.0∕(bar)

wX-q

=1.0∕(1000 m

∕h).

F.Hennings et al.

1 3

solutions for the problem. This algorithm also takes care of the restrictions regarding

the operation mode transition times from Sect.2.8, which were not included in

5.2 Computational setup

We performed our computations on a cluster using 4 cores and 16GB of RAM

of a machine composed of two Intel Xeon CPU E5-2670 v2 running at 2.50GHz.

As a solver for the underlying MIP problems, we used Gurobi in version 8.1.0,

see GurobiOptimization (2018), which we accessed via the Pyomo modeling lan-

guage introduced by Hart etal. (2011, (2017). Since the corresponding MIP models

were numerically challenging, we used the solver with the maximal NumericFocus

parameter. Furthermore, we specified the optimality conditions and maximum run

times for each solved model variant as

Finally, we chose the rolling horizon parameter h to be 4, i.e., we solved the tran-

sient smoothing for 4 future time steps. We found in our experiments that this num-

ber represents a good trade-off between efficient model solving speed and foresight-

edness of the obtained solution.

5.3 Results

As described in Sect. 5.1 above, we tested Algorithm 1 on

159 start times×

7 network stations ×4 time horizons =4452

instances. For all of these instances, it

found a feasible solution. In particular, we always found a feasible solution via the

initial solution creation in Algorithm2 described in Sect.4.2.1, although the algo-

rithm itself does not guarantee this.

In Fig.5, one can find an overview of the average run times of Algorithm1 sorted

by network station and number of solved time steps. As one can see, the ranking

of the single stations regarding run time as well as the ratio between the single sta-

tions run times is stable over the different time horizon partitions. Furthermore, the

run time per station corresponds roughly to the complexity deduced from the station

statistics in the sense that more nodes, arcs, and operation modes make the overall

problem more complex to solve and therefore increase the run time. The exception

to this reasoning is station C, which has a rather high average run time, although its

node and arc counts are comparable to station D. Its number of operation modes is

even the smallest of all seven stations. We assume that this is due to the station’s

topology, in which the regulators may be arranged in a way that is hard to solve for

the final smoothing step of the algorithm. The higher portion of runtime spent in the

smoothing of station C in comparison to station D is illustrated in Fig.7. We have

been unable to prove this assumption so far, however. Regarding model complex-

ity, we did not observe that the existence of a non-empty set

of flow direction

conditions adds any notable difficulty to solving station E, which is the only station

Variant Rel. Gap Abs. Gap Time limit

s,Psf 1E-4 1E-2 10

5E-3 1E-2 60s.

1 3

Controlling transient gas flow inreal‑world pipeline…

Table 2 The statistics of model

for each network station and each number of time steps. The triples state the number of continuous variables, binary variables, and con-

straints, each rounded to the nearest hundred. Note that

does not include the restrictions regarding the operation mode transition times from Sect.2.8

Station Number of time steps

12 24 48 96

Cont Binary Constr Cont Binary Constr Cont Binary Constr Cont Binary Constr

A 2200 900 7300 4300 1700 14400 8500 3400 28,600 17,000 6800 57,100

B 2000 700 6900 4000 1300 13700 8000 2600 27,300 15,900 5200 54,400

C 3100 700 5000 6100 1500 10000 12,100 2900 19900 24,200 5900 39,800

D 2500 1700 8100 4900 3400 16000 9700 6800 31,800 19,300 13,600 63,500

E 5100 2000 11,800 10,100 3900 23,500 20,200 7800 46,800 40,300 15,700 93,500

F 6000 19700 117,300 11,900 36,600 231,600 23,700 70,400 460,100 47,300 137,900 917,500

G 10,200 18600 67,400 20,300 35,900 132,600 40,300 70,500 265,200 80,500 139,700 528,800

F.Hennings et al.

1 3

incorporating these constraints. The last fact we want to highlight here is the almost

perfectly proportional average run time increase with increasing number of time

steps. As mentioned in Sect.4.3, this property was an explicit goal when designing

the algorithm.

In addition to the average run times, we also depict the distribution of run times

sorted by station in Fig.6 for instances solved for 12 time steps. The distribution of

the other time steps are quite similar in appearance and can be found in the Appen-

dix as Figs.9, 10, and 11.

For the fast solving stations A, B, and D, the run times of the instances are all

very similar in that the lines representing the 25th percentile and the 75th percentile

are extremely close. Station C is again outstanding by having a rather large spread

of run times only matched by the biggest and most complex station G. However, in

general the maximum run time is less than an order of magnitude larger than the

median run time for all of the stations. This means that for our instance set, even the

extreme cases are still in reach and have a somewhat similar run time to the rest of

the instances, making Algorithm1 well suited for strict time limit situations, which

we would face in a production environment.

As the last graphic related to the run time, Fig.7 shows the average portion of

run time spent in each individual subroutine of Algorithm1 for 12 time steps as well

as for 96 time steps. It states that the transient smoothing dominates the run time

already for 12 time steps. This run time domination increases as the number of time

steps do, while the heuristic’s run time proportion, and therefore impact, decreases.

This can also be verified through the 24 and 48 time step instances displayed in

Fig.12 in the Appendix. The observed time distributions fit the structure of the sin-

gle subroutines from Sect.4. While in the initial solve, as well as in the smooth-

ing, more instances of the harder model variants

and

are solved with each

additional time step, the improvement heuristic iterates over pairs of subsequent

time steps with different operation modes. Since the differently sized time horizon

Fig. 5 Average run times for Algorithm1 sorted by number of time steps and network station

1 3

Controlling transient gas flow inreal‑world pipeline…

partitions all cover the same 12 hours, the number of total operation mode changes

in the result of the initial solve is very unlikely to increase much with an increasing

time granularity. This is also due to the initial solve algorithm itself, which tries to

use the same operation mode as long as possible before changing to a different one.

Apart from the run time, we were also interested in the quality of the solution of

Algorithm1. Therefore, we solved the complete model

directly on the smaller

network stations A to E to obtain a valid lower bound. Note that

itself is already a

Fig. 6 Run times of instances with 12 time steps, sorted by station. One dot represents one instance,

however, dots may overlap when too many have similar run time. In addition, we draw a typical whisker-

less box plot. Therefore the lines represent the 25th percentile, the median, and the 75th percentile

Fig. 7 Portion of run time spent in the single subroutines of Algorithm1 displayed for each of the single

stations. Figure7a represents the results for 12 time steps, Fig.7b the results for 96 time steps

F.Hennings et al.

1 3

relaxation of the problem solved by Algorithm1 since the transition time restric-

tions for operation modes introduced in Sect.2.8 are not included in this model. The

obtained lower bound can be used to give an upper bound on the relative difference

of the solution found by Algorithm1 to the optimal solution. We call this difference

gap and define it as

obj(solution)−lowerBound

obj(solution)

. Since the lower bound cannot be negative,

the gap is always smaller than or equal to 100%.

We solved model

on stations A to E for 12 time steps using the same computa-

tional setup as described above with a time limit of 10 hours. Furthermore, we gave the

solution obtained from Algorithm1 as a starting solution. The solver finished for all

159 instances of the stations A and B within the time limit. For the stations C, D, and

E, we hit the time limit for 47, 8, and 18 instances, respectively, resulting in potentially

sub-optimal lower bounds in these cases. For the stations F and G, more than half of the

instances did not finish in time, which is why we did not include these in our evaluation.

The results can be found in Fig.8, where we plotted the gap for all 5 network

stations on a logarithmic scale and then plotted the same data for stations C and E

Fig. 8 Gap between the solution of Algorithm1 and a lower bound obtained from solving

directly.

We considered the stations A to E and solved them using 12 time steps. One dot represents one instance,

however, dots may overlap when too many have similar run time. The instances with potentially subop-

timal lower bounds due to hitting the time limit are marked by triangles instead of dots. In addition, we

draw a typical whiskerless box plot. Therefore the lines represent the 25th percentile, the median, and the

75th percentile. The left graphic uses a logarithmic scale for the gap, and the right one displays the same

data for stations C and E on a linear scale. For all instances in which the objective function value and

the lower bound are smaller than 0.1, we defined the gap to be zero. Furthermore, we plotted all values

below

10−3

for the logarithmic scale

1 3

Controlling transient gas flow inreal‑world pipeline…

on a linear scale to better depict their distributions. Note that the instances that have

a potentially suboptimal lower bound are marked with triangles instead of dots. As

one can observe, we have good results for all of the network stations, with more than

half of the instances for each station having a solution with at most 10% difference

to the lower bound. For the best three stations A, B, and D, we can even state that

at least 75% of the values have less than 1% difference to the optimal solution. The

picture is more diverse for the other two stations, in which the 75% percentile is

between 20% and 30% gap. However, for those stations, we have the highest number

of instances for which

did not solve the problem in time. Therefore, there is the

potential for further improvement of the gap by increasing the lower bound here.

Most of the instances with potentially suboptimal lower bounds also have high gap

values. The picture is extreme for Station C, where all but 3 instances above the 75%

percentile had suboptimal lower bounds, making it very likely that one can further

decrease the gap value by increasing the run time limit of

directly.

Altogether, the results show that Algorithm1 is able to find good solutions reli-

ably and fast. Hence it can solve even large instances of the presented gas transpor-

tation problem, which would be far out of reach when trying to solve the original

MIP model formulation

directly.

A comparison to experiments done in the existing literature on transient gas net-

work optimization confirms our claim of a competitive algorithm, albeit the fact that

a detailed evaluation of the different approaches is not possible since every publica-

tion considers a different model. However, we are in general able to solve our tran-

sient gas transportation problem faster than the algorithms presented in the literature

solve their corresponding problem versions. This is true, even though we consider

larger instances. All the experiments conducted by Domschke etal. (2011), Hahn

etal. (2017), and Burlacu etal. (2019) consider networks of at most 20 nodes and

arcs. Furthermore, the observed run times and given time limits are often larger than

those of our algorithm. For the largest instance in Domschke etal. (2011) solved

using 4 time steps, their algorithm hits the time limit of 10000 seconds with a gap of

34%. In Hahn etal. (2017), the time horizon was split into 180 time steps, while in

only 10 of these time steps discrete control decisions were allowed. There, the solv-

ing times ranged from 2 hours up to the time limit of 50 hours. Finally, the largest

instance considered in Burlacu etal. (2019) used 48 time steps, and their algorithm

hits the time limit of 3 hours with a gap of 5.5%. However, as already mentioned

above, the problems described in these publications differ from the one presented in

this paper. In general, each of them considers a more accurate, nonlinear pipe model.

At the same time, they only use simple models for the compressors consisting of sin-

gle compressor units without configurations or interdependent modes. Furthermore,

we included a series of realistic constraint types, like the choice of flow directions

together with the corresponding exit pressure constraints and flow direction condi-

tions, or the operation mode transition times. For these, only Burlacu etal. (2019)

included a similar type of constraint by introducing switching times for valves and

compressor units.

F.Hennings et al.

1 3

6 Conclusion

We presented the transient gas network transportation problem on network stations,

which represent the intersection points of major transportation pipelines and con-

tain the majority of active elements to control the network. For this problem, we

introduced a mixed-integer programming model, including a sophisticated model for

compressor stations, as well as additional variables and constraints for the network

station itself. For the pipes, we found that due to their shortness, they have less over-

all impact in network stations, which enabled us to use a linear description for them.

The MIP model was observed, through the use of state-of-the-art solvers, to be

intractable for large network stations or a large number of time steps. Therefore, we

developed a specialized algorithm to solve the problem. Its design is influenced by

the fact that network stations contain only short pipes and hence a negligible amount

of linepack. This causes the transient problem to be a “stationary-flavored” one. For

this reason, our algorithm determines the operation modes and gas flow directions

of the network station based on solving multiple stationary (or steady-state) vari-

ants of the MIP. By using this approach, we were also able to satisfy the transition

time restrictions for the operation modes, which we excluded from the original MIP

formulation. In order to obtain a feasible solution for the problem, we finally solved

another variant of the MIP with fixed operation modes and flow directions for the

network station in a rolling horizon fashion.

To verify the competitiveness of our algorithm regarding run time, reliability to

find a solution, and quality of that solution, we performed tests on 159 different sce-

narios based on past flow situations in 7 real-world network stations provided by our

project partner. Our algorithm was able to compute feasible solutions for all of the

presented instances very fast. Even for the largest of the presented network stations

consisting of more than 100 nodes and arcs as well as over 1000 different network

station operation modes, our algorithm terminates on average in under 20 minutes

for each of the stations. By running experiments using 4 different types of granular-

ity, we found that the run time increases proportionally to the number of time steps

used, indicating that our algorithm scales very well. In terms of quality of the solu-

tion, we tested the results against a lower bound obtained from solving the original

MIP formulation for 12 time steps on all except the largest two network stations. For

the worst of the five remaining stations, more than half of the instances had a solu-

tion with at most 10% difference to the lower bound. For the best three stations, we

found near-optimal solutions with less than 1% difference for more than 75% of the

instances. Altogether we were able to show that our algorithm can reliably find good

solutions to the problem in a short amount of time.

There are a lot of different possibilities to continue this research. To increase the

model accuracy, the approximative linearization of the friction in pipes and resis-

tors, as well as the maximum power bound constraint of compressor units could

be replaced by their original nonlinear versions. This would turn the model into

a MINLP and thereby increase its complexity. From a theoretical point of view,

extending Algorithm 2 such that it has the theoretical guarantee to always find

1 3

Controlling transient gas flow inreal‑world pipeline…

existing feasible solutions would improve its robustness. Finally, we name ramp-

up and cool-down times for compressor units, as well as using target value based

control of regulators and compressor units, as examples of extensions to our model.

These are examples of a never-ending list of unique elements and extra constraints

that shapes the complicated business of real-world gas network operation.

Acknowledgements The work for this article has been conducted in the Research Campus MODAL

funded by the German Federal Ministry of Education and Research (BMBF) (fund number 05M14ZAM)

Funding Open Access funding enabled and organized by Projekt DEAL.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License,

which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as

you give appropriate credit to the original author(s) and the source, provide a link to the Creative Com-

mons licence, and indicate if changes were made. The images or other third party material in this article

are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the

material. If material is not included in the article’s Creative Commons licence and your intended use is

not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission

directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen

ses/by/4.0/.

Appendix

See Table3; Figs.9, 10, 11 and 12.

F.Hennings et al.

1 3

Table 3 List of all used variables, specifying their domain, meaning and unit. All variables are defined

for

t∈T

, not taking into account the given values for the initial state, see Sect. 2.9. Note that

1 bar =105Pa

Variable Meaning Unit

pv,t

∈ℝ≥0

Pressure at node

v∈V

bar

qa,t

∈ℝ

Flow over arc

a∈

⧵

kg/s

qv,a,t

∈ℝ

Flow into or out of pipe

a∈Api

kg/s

dv,t

∈ℝ

Inflow into boundary node

v∈Vb

kg/s

mcf

c,a,t

∈{0, 1}

Selection of configuration

c∈Ca

for

a∈Acs

mby

a,t

∈{0, 1}

Selection of bypass mode for

a∈Arg ∪Acs

mcl

a,t

∈{0, 1}

Selection of closed mode for

a∈Arg ∪Acs

mac

a,t

∈{0, 1}

Selection of active mode for

a∈Arg

mop

a,t

∈{0, 1}

Selection of open/closed mode for

a∈Ava

pby

a,t

∈ℝ≥0

Pressure in bypass mode for

a∈Acs

bar

qby

a,t

∈ℝ

Flow in bypass mode for

a∈Acs

kg/s

pl-cl

a,t

∈ℝ≥0

Incoming pressure in closed mode for

a∈Acs

bar

pr-cl

a,t

∈ℝ≥0

Outgoing pressure in closed mode for

a∈Acs

bar

l-cf

c,a,t

∈ℝ≥0

Incoming pressure in configuration

c∈Ca

a∈Acs

bar

pr-cf

c,a,t

∈ℝ≥0

Outgoing pressure in configuration

c∈Ca

a∈Acs

bar

qcf

c,a,t

∈ℝ

Flow in configuration

∈

a∈Acs

kg/s

mom

o,t

∈{0, 1}

Selection of mode

o∈O

mfd

f,t

∈{0, 1}

Selection of flow direction

f∈F

𝜎p+

v,t

∈ℝ≥0

Positive pressure slack for boundary node

v∈Vb

bar

𝜎p−

v,t

∈ℝ≥0

Negative pressure slack for boundary node

v∈Vb

bar

𝜎d+

v,t

∈ℝ≥0

Positive flow slack for boundary node

v∈Vb

kg/s

𝜎d−

v,t

∈ℝ≥0

Negative flow slack for boundary node

v∈Vb

kg/s

𝛿om

∈{0, 1}

Operation mode change 1

𝛿rg

a,t

∈{0, 1}

Mode change for regulator

a∈Arg

𝛿us

u,t

∈{0, 1}

Start of compressor unit

u∈Ua

for

a∈Acs

𝛿rg-pl

a,t

∈ℝ≥0

Change of incoming pressure of active

a∈Arg

bar

𝛿rg-pr

a,t

∈ℝ≥0

Change of outgoing pressure of active

a∈Arg

bar

𝛿rg-q

a,t

∈ℝ≥0

Change of flow over of active

a∈Arg

kg/s

𝛿cs-pl

a,t

∈ℝ≥0

Change of incoming pressure of active

a∈Acs

bar

𝛿cs-pr

a,t

∈ℝ≥0

Change of outgoing pressure of active

a∈Acs

bar

𝛿cs-q

a,t

∈ℝ≥0

Change of flow over of active

a∈Acs

kg/s

1 3

Controlling transient gas flow inreal‑world pipeline…

Fig. 9 Run times of instances with 24 time steps, sorted by station. One dot represents one instance,

however, dots may overlap when too many have similar run time. In addition, we draw a typical whisker-

less box plot. Therefore the lines represent the 25th percentile, the median, and the 75th percentile

Fig. 10 Run times of instances with 48 time steps, sorted by station. One dot represents one instance,

however, dots may overlap when too many have similar run time. In addition, we draw a typical whisker-

less box plot. Therefore the lines represent the 25th percentile, the median, and the 75th percentile

F.Hennings et al.

1 3

References

Balas E (1985) Disjunctive programming and a hierarchy of relaxations for discrete optimization prob-

lems. SIAM J Algeb Discrete Methods 6(3):466–486. https ://doi.org/10.1137/06060 47

Balas E (2018) The convex hull of a disjunctive set. Disjunctive programming. Springer, Cham, pp 17–39

Bonami P, Lodi A, Tramontani A, Wiese S (2015) On mathematical programming with indicator con-

straints. Math Prog 151(1):191–223. https ://doi.org/10.1007/s1010 7-015-0891-4

Brouwer J, Gasser I, Herty M (2011) Gas pipeline models revisited: model hierarchies, nonisother-

mal models, and simulations of networks. Multisc Model Simul 9(2):601–623. https ://doi.

org/10.1137/10081 3580

Fig. 11 Run times of instances with 96 time steps, sorted by station. One dot represents one instance,

however, dots may overlap when too many have similar run time. In addition, we draw a typical whisker-

less box plot. Therefore the lines represent the 25th percentile, the median, and the 75th percentile

Fig. 12 Portion of run time spent in the single subroutines of Algorithm1 displayed for each of the sin-

gle stations. Figure12a represents the results for 24 time steps, Fig.12b the results for 48 time steps

1 3

Controlling transient gas flow inreal‑world pipeline…

Burlacu R, Egger H, Groß M, Martin A, Pfetsch M, Schewe L, Sirvent M, Skutella M (2019) Maximizing

the storage capacity of gas networks: a global MINLP approach. Optim Eng 20(2):543–573. https ://

doi.org/10.1007/s1108 1-018-9414-5

Domschke P, Geißler B, Kolb O, Lang J, Martin A, Morsi A (2011) Combination of nonlinear and linear

optimization of transient gas networks. INFORMS J Comput 23(4):605–617

Domschke P, Hiller B, Lang J, Tischendorf C (2017) Modellierung von Gasnetzwerken: Eine Übersicht

Ehrhardt K, Steinbach MC (2003) Nonlinear optimization in gas networks. Technical Report 03-46, ZIB,

Takustr. 7, 14195 Berlin

Federal Ministry for Economic Affairs and Energy (2019) Still indispensable for a reliable energy supply.

https ://www.bmwi.de/Redak tion/EN/Dossi er/conve ntion al-energ y-sourc es.html. Accessed 31 Oct

2019

Fügenschuh A, Geißler B, Gollmer R, Morsi A, Pfetsch ME, Rövekamp J, Schmidt M, Spreckelsen K,

Steinbach MC (2015) Physical and technical fundamentals of gas networks. In: Koch etal. (2015)

Gugat M, Leugering G, Martin A, Schmidt M, Sirvent M, Wintergerst D (2018) MIP-based instantaneous

control of mixed-integer PDE-constrained gas transport problems. Comput Optim Appl 70(1):267–

294. https ://doi.org/10.1007/s1058 9-017-9970-1

GurobiOptimization L (2018) Gurobi optimizer reference manual, version 8.1.0. http://www.gurob i.com

Hahn M, Leyffer S, Zavala VM (2017) Mixed-integer PDE-constrained optimal control of gas networks.

Preprint

Hart WE, Watson JP, Woodruff DL (2011) Pyomo: modeling and solving mathematical programs in

Python. Math Prog Comput 3(3):219–260

Hart WE, Laird CD, Watson JP, Woodruff DL, Hackebeil GA, Nicholson BL, Siirola JD (2017) Pyomo-

optimization modeling in python, vol 67, 2nd edn. Springer, Berlin

Hennings F (2018) Benefits and limitations of simplified transient gas flow formulations. In: Operations

research proceedings 2017. Springer, pp 231–237

Hoppmann K, Hennings F, Lenz R, Koch T (2019) Optimal operation of transient gas transport networks.

Technical Report 19-23, ZIB, Takustr. 7, 14195 Berlin

Humpola J, Fügenschuh A, Hiller B, Koch T, Lehmann T, Lenz R, Schwarz R, Schweiger J (2015) The

specialized MINLP approach. In: Koch etal. (2015)

Koch T, Hiller B, Pfetsch ME, Schewe L (eds) (2015) Evaluating gas network capacities, MOS-SIAM

series on optimization, vol21. SIAM

Kolb O, Lang J, Bales P (2010) An implicit box scheme for subsonic compressible flow with dissipative

source term. Numeri Algorith 53(2–3):293–307

Mahlke D, Martin A, Moritz S (2007) A simulated annealing algorithm for transient optimization in gas

networks. Math Methods Oper Res 66(1):99–115. https ://doi.org/10.1007/s0018 6-006-0142-9

Mak TW, VanHentenryck P, Zlotnik A, Hijazi H, Bent R (2016) Efficient dynamic compressor optimiza-

tion in natural gas transmission systems. In: American control conference (ACC), 2016, IEEE, pp

7484–7491. https ://doi.org/10.1109/ACC.2016.75268 55

Moritz S (2007) A mixed integer approach for the transient case of gas network optimization. PhD thesis,

Technische Universität Darmstadt, Darmstadt

Nikuradse J (1950) Laws of flow in rough pipes. National Advisory Committee for Aeronautics

Washington

Odom FM, Muster GL (2009) Tutorial on modeling of gas turbine driven centrifugal compressors. In:

PSIG-09A4, Pipeline Simulation Interest Group

Osiadacz AJ (1996) Different transient flow models-limitations, advantages, and disadvantages. In: PSIG-

9606, Pipeline Simulation Interest Group

Papay (1968) A termeléstechnológiai paraméterek változása a gáztelepek muvelése során. OGIL Musz

Tud Kozl

Pfetsch ME, Fügenschuh A, Geißler B, Geißler N, Gollmer R, Hiller B, Humpola J, Koch T, Lehmann

T, Martin A, Morsi A, Rövekamp J, Schewe L, Schmidt M, Schultz R, Schwarz R, Schweiger J,

Stangl C, Steinbach MC, Vigerske S, Willert BM (2014) Validation of nominations in gas network

optimization: models, methods, and solutions. Optim Methods Softw. https ://doi.org/10.1080/10556

788.2014.88842 6

Ríos-Mercado RZ, Borraz-Sánchez C (2015) Optimization problems in natural gas transportation sys-

tems: a state-of-the-art review. Appl Energy 147:536–555

Saleh JM (2002) Fluid flow handbook. McGraw-Hill, London

Walther T, Hiller B (2017) Modelling compressor stations in gas networks. Technical Report 17-67, ZIB,

Takustr. 7, 14195 Berlin

F.Hennings et al.

1 3

Zlotnik A, Chertkov M, Backhaus S (2015) Optimal control of transient flow in natural gas networks. In:

54th IEEE conference on decision and control (CDC), IEEE, pp 4563–4570. https ://doi.org/10.1109/

CDC.2015.74029 32

Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published

maps and institutional affiliations.

Affiliations

FelixHennings1,2 · LovisAnderson1 · KaiHoppmann‑Baum1,2 ·

MarkTurner1,2 · ThorstenKoch2,1

Lovis Anderson

Kai Hoppmann-Baum

Mark Turner

Thorsten Koch

1 Zuse Institute Berlin, Takustr. 7, 14195Berlin, Germany

2 Chair ofSoftware andAlgorithms forDiscrete Optimization, Institute ofMathematics,

Technische Universität Berlin, Straße des 17. Juni 135, 10623Berlin, Germany