Abstract
In this paper, we describe an algorithmic framework for the optimal operation of transient gas transport networks consisting of a hierarchical MILP formulation together with a sequential linear programming inspired post-processing routine. Its implementation is part of the KOMPASS decision support system, which is currently used in an industrial setting. Real-world gas transport networks are controlled by operating complex pipeline intersection areas, which comprise multiple compressor units, regulators, and valves. In the following, we introduce the concept of network stations to model them. Thereby, we represent the technical capabilities of a station by hand-tailored artificial arcs and add them to network. Furthermore, we choose from a predefined set of flow directions for each network station and time step, which determines where the gas enters and leaves the station. Additionally, we have to select a supported simple state, which consists of two subsets of artificial arcs: Arcs that must and arcs that cannot be used. The goal is to determine a stable control of the network satisfying all supplies and demands. The pipeline intersections, that are represented by the network stations, were initially built centuries ago. Subsequently, due to updates, changes, and extensions, they evolved into highly complex and involved topologies. To extract their basic properties and to model them using computer-readable and optimizable descriptions took several years of effort. To support the dispatchers in controlling the network, we need to compute a continuously updated list of recommended measures. Our motivation for the model presented here is to make fast decisions on important transient global control parameters, i.e., how to route the flow and where to compress the gas. Detailed continuous and discrete technical control measures realizing them, which take all hardware details into account, are determined in a subsequent step. In this paper, we present computational results from the KOMPASS project using detailed real-world data.
Similar content being viewed by others
1 Introduction
Natural gas is and will remain one of the major energy sources in Europe (https://www.bmwi.de/Redaktion/EN/Dossier/conventional-energy-sources.html). Furthermore, it is often considered an important transit medium towards a low- or no-carbon future (Winegarden 2019). While the overall gas consumption in Germany is assumed to remain constant in the future, the hourly supplies and demands at the sources and sinks of the network are expected to become more volatile. An example reason for this behaviour is the growing usage of renewable energy, e.g., solar and wind power. While their share in the energy mix is going to increase due to the planned nuclear and coal phase-outs, there is a lot of uncertainty regarding their production. One possibility to mitigate this effect is to use natural gas fired power plants, which can be ramped up on short notice. Hence, for the Open Grid Europe GmbH (OGE) (https://www.open-grid-europe.com/), one of the largest transport system operators (TSOs) in Germany (https://www.open-grid-europe.com/cps/rde/oge-internet/hs.xsl/Strukturdaten-gemass-27-Abs-2-GasNEV-654.htm?rdeLocaleAttr=de), a more robust, secure, and stable control of the network becomes inevitable in order to guarantee security of supply. Thus, the idea for KOMPASS (Kontinuierliches Optimierungsmodul zur prognoseabgesicherten Systemsteuerung/Continuous optimization module for a prognosis-based system control), a decision support system for the dispatchers, who control the gas network, was born and realized within the GasLab of the Research Campus MODAL (http://forschungscampus-modal.de/). Before we explain the architecture of KOMPASS in detail in Sect. 2 and thereby motivate the model presented in this paper, we first give an overview of important related work.
The optimization of natural gas transport through pipeline networks is a challenging task due to two crucial aspects: The physics of gas flow and the combinatorics behind the setup of compressor units together with imposed technical restrictions and limitations. The gas flow through cylindric pipelines is commonly described by the so-called Euler equations (Osiadacz 1996), a set of non-linear hyperbolic partial differential equations (PDEs), see Sect. 4.8 for details. On the other hand, compressor units can run sequentially or in parallel by opening and closing surrounding valves in order to achieve the required compression ratio and flow rate. Additionally, they feature feasible operating ranges and are subject to a non-linear power bound. More detailed explanations on these topics can be found in Hennings et al. (2019) and Koch et al. (2015).
While the stationary case, i.e., determining a feasible network state given the necessary boundary values, has gained a lot of attention in recent years, see Koch et al. (2015) for an extensive overview, research regarding the transient case is still in the early stages. One of the first papers on transient gas transport optimization to appear is Moritz (2007). Here, a mixed integer linear program (MILP) featuring independent single compressor units is introduced. Furthermore, the gas flow in pipelines is modelled using piecewise linear functions. Pure non-linear programs (NLPs) are considered in Mak et al. (2016) and Zlotnik et al. (2015), which decide on the compression ratios of compressors, while minimizing fuel consumption. Recently Gugat et al. (2018) and Burlacu et al. (2019) have been published, which both make use of special discretization schemes for the Euler equations and again consider independent single compressor units. Gugat et al. (2018) use a linear feasible region for the compressors and minimize the deviation from future flow and pressure values by iteratively solving a MILP and an NLP model for each time step. In contrast, Burlacu et al. (2019) impose lower and upper bounds on the compression ratios and on the achievable pressure differences of the compressor units. Here, the amount of gas stored in the network is maximized by alternatingly solving a MILP and an NLP model.
On the other hand, due to the liberalization of many energy markets, where so-called entry-exit models were introduced, many new optimization problems arose in these areas. For such problems, hierarchical optimization turns out to be well suited for modelling purposes, and this is not only because, for example, the entry-exit model for the European natural gas market (Rövekamp 2015) itself can be described using a multi-level formulation (Grimm et al. 2019). For the US market a discrete bi-level programming approach is used to solve a cash-out problem, where a gas shipper weighs daily delivery imbalances against penalties claimed by the transport company (Dempe et al. 2005; Kalashnikov and Ríos-Mercado 2002). Bi- and multi-level programming in general has recently experienced a renaissance and is used more often to model a variety of real-world applications. Areas include but are not limited to network design (Gao et al. 2005), capacity planning (Garcia-Herreros et al. 2016), toll setting (Brotcorne et al. 2001), robust unit commitment (Jiang et al. 2011), or critical infrastructure protection (Alderson et al. 2011). Detailed overviews on possible applications can be found in Colson et al. (2007), Kalashnikov et al. (2015), and Lu et al. (2016), and for an introduction to hierarchical and bi-level optimization we refer to Dempe (2002) and Migdalas et al. (2013).
The remainder of this paper is structured as follows: In Sect. 2 we give an overview of the architecture of the KOMPASS decision support system to further motivate the problem formulation presented in this paper. Afterwards, we introduce the concept of network stations, a modelling approach to simplify complex pipeline intersection areas, and explain the corresponding mathematical formulation in Sect. 3. Subsequently, we introduce our tri-level transient gas flow MILP model in Sect. 4, for which a solution approach is presented in Sect. 5. Afterwards, we introduce a sequential linear programming inspired post-processing routine to derive physically more accurate solutions in Sect. 6. The solution approach and the post-processing routine from the last two sections are subsequently combined in an algorithmic framework, which is presented in Sect. 7. We conclude with computational experiments regarding this framework, which were conducted on real-world instances, and an outlook on future improvements and extensions in Sects. 8 and 9, respectively.
2 KOMPASS—A decision support system for gas network control
The dispatchers at OGE control the gas network mainly based on their personal set of skills, e.g., knowledge from training they receive and experience. However, since they are starting to face more volatile and previously unseen transport scenarios, guaranteeing a safe operation and security of supply has become more difficult. Hence, the idea for a decision support system for the control of transient gas transport networks was born: KOMPASS (Kontinuierliches Optimierungsmodul zur prognoseabgesicherten Systemsteuerung/Continuous optimization module for a prognosis-based system control). Its architecture, as currently implemented and running at OGE, is shown in Fig. 1.
Architecture of KOMPASS and flow of information within it. Red squares denote computational modules, blue parallelograms describe input data (if in top row) or intermediate output data (if not in top row), which serve as input for subsequent algorithmic modules, and green ovals denote the output data
First of all, KOMPASS receives the topology of the network, its current state, and prognosis data as input. The latter features historic gas flows at the sources and sinks of the network, weather data, and information regarding work- and holidays. Based on this, supplies and demands are predicted using the approach presented in Petkovic et al. (2019), and time series of future pressure values for the sources are subsequently determined by a heuristic.
As output, highly detailed technical control recommendations shall be determined for the next 12 h and all remotely controllable elements. To be able to create such recommendations, we are given additional information describing these elements. Thus, their characterizations are part of the input and contained in the intersection area data, which exists for each complex pipeline intersection area. Another important part of the intersection area data are the operation modes. Each operation mode features detailed technical settings for all remotely controllable elements at the corresponding junction and thereby enables certain technical control possibilities. A main goal of the KOMPASS decision support system is to determine operation modes for each junction such that the given supplies and demands are satisfied while the stability of the network is maximized. Switching from one operation mode to another is called a technical measure and we consider the operation of a network to be more stable the less technical measures have to be applied in the following, see Hennings et al. (2019) for details.
In the real-world the dispatchers try to control the network using technical measures only. But since this is not always possible, they additionally have some non-technical measures at hand. The most common and standardized ones result in changes to the supplies and demands by for example either buying or selling gas, i.e., the usage of so-called balancing energy, or by using contractual options like the interruption of customers. If changing the supplies and demands alone does not guarantee a feasible control of the network, the last option is to ask neighboring TSOs for future pressure changes at some sources. In practice, this is done by phone calls and it can therefore be seen as an emergency and non-standardized option. This establishes a hierarchy on the usage of measures. If the network can be controlled by technical measures only, i.e., without the usage of non-technical measures, this is most favorable. As a second option, deviations from the supplies and demands are allowed. If there still does not exist any feasible control for the network, the last option consists of a change of future source pressures.
Thus, given the network topology, its current state, time series on supplies and demands as well as on the pressures at the sources, and the intersection area data, we have to solve a transient gas network control problem on a large-scale and complex real-world gas network within KOMPASS. This problem formulation has to first of all incorporate the operation modes, which include the complex combinatorics of the setup of compressor stations. Second, an as accurate as possible physical model of the transient gas flow in pipelines. And third, the described hierarchy of technical and non-technical measures. Additionally, due to the nature of the application for which KOMPASS is designed, its run time does play an important role. The decision support system is supposed to run 24/7 and continuously provide technical and non-technical control recommendations to the dispatchers.
However, our first experiments showed that models incorporating all these needs simultaneously were computationally intractable or not solvable within any reasonable amount of time. Thus, we decided to split the complexity and introduced a two-stage approach. First, a transient control problem using hand-tailored simplified models for each of the complex pipeline intersection areas is solved. The original junctions are replaced by these simplified graph representations, which we call network stations in the following, and we call the resulting network macroscopic. Their description is also part of the intersection area data and their derivation is explained in more detail in Sect. 3. Furthermore, non-technical measures are recommended based on the results of this model.
Afterwards, the resulting flow and pressure values at the boundaries of the network stations serve as input for highly detailed models on the original complex pipeline intersection areas. In particular, these models validate whether there exist actual operation modes realizing the given pressure and flow scenarios (Hennings et al. 2019). Here, stationary models focusing on the combinatorics and the technical restrictions of compressor stations are solved in a first step. The rationale behind this is that intersection areas contain only pipelines of short length, which cannot store or provide much gas for future usage, i.e., they do not feature a lot of linepack. Therefore the transient aspect of the gas flow can be neglected. However, this aspect is included in a second step, where a corresponding mathematical model is solved using a rolling horizon approach. For more details on the second stage of KOMPASS we refer to the original source (Hennings et al. 2019).
The paper at hand deals with the first stage of the approach, i.e., formulating and solving the transient gas flow problem on macroscopic networks. Therefore, we introduce a tri-level MILP formulation to accommodate for the hierarchy regarding the technical and non-technical measures. A top to bottom explanation of the rationale behind our hierarchical modelling approach is depicted in Fig. 2.
Top to bottom explanation of the rationale behind our tri-level formulation. Note that change instructions are executed such that the sums of the absolute deviations from the given input values are minimized
However, for the sake of comprehensibility, we additionally give a bottom to top explanation in the following. The third level features the technical control problem, which tries to maximize the stability of the network by minimizing changes in the operation of the network stations, i.e., by minimizing the usage of technical measures. The second level minimizes the changes of supplies and demands w.r.t. the sum of absolute deviations that is necessary to guarantee feasibility of the third level. The first level pursues a similar goal but minimizes the sum of deviations from future pressure values at the sources instead. Thereby, it takes all possible actions of the second level into account. In other words, the first and the second level minimize the extent of their corresponding non-technical measures in hierarchical order to ensure that a feasible technical control of the network exists.
3 Network station fundamentals
A majority of the remotely controllable elements within gas transport networks, i.e., elements whose behaviour can be remotely controlled by the dispatchers, such as compressor stations, regulators, and valves, is located at intersections of major transportation pipelines. For each of these junctions and each point in time, exactly one operation mode is in use. Due to the amount of possible operation modes and the induced complexity, together with experts from OGE we developed the hand-tailored simplified graph representations called network stations in order to summarize and approximate the technical control capabilities. While detailed mathematical formulations and explanations regarding network stations can be found in Sect. 4.9, we briefly explain the basic idea and the process of their derivation here. To do this, an example of a simplification is shown in Fig. 3.
The colored triangles represent sources (green triangle) and sinks (red triangle), which are closely located to the network station. Further, the other single network elements are depicted as 
(pipe), 
(valve/shortcut), 
(regulator/regulating arc), 
(compressor station/compressing arc), and 
(bi-directed regulating arc). The fence groups of the network station are highlighted by colored circles
As mentioned before, the derivation of the network station models is currently manually done by experts at OGE, who know their network, its elements and its control very well. Thus, we can only give some intuition of how they are created here by using examples. But as we discuss later in the outlook in Sect. 9, we are currently working on an automatized process to derive them.
First of all, the intersection areas are identified as connected subgraphs of the network. Their layouts are created with the goal in mind to include as many remotely controllable elements as possible while only containing few pipelines of significant length. The latter follows from the idea behind the two-stage approach in KOMPASS, where the transient behaviour of the gas flow is considered less important in the second stage, see Sect. 2 for the details.
The nodes at the boundaries of these subgraphs are called fence nodes. If a subset of the fence nodes features the same behaviour, e.g., all are connected to pipelines of large diameter, which run in parallel and nearly always possess the same pressure level as well as the same direction and amount of flow, they are merged in the network station topology and called a fence group. In Fig. 3 the fence groups are indicated using colored circles.
In the next step we create the topology of the network station by removing the interior of the subgraph and adding auxiliary nodes together with artificial arcs, which connect them and the fence nodes. There are four types of artificial arcs: Shortcuts, which can be seen as the equivalents of valves, regulating arcs, which can be seen as regulators, compressor arcs, which shall capture the pressure increasing capabilities of compressor stations, and combined arcs, which can work as either regulating or compressor arcs. Further, for each of these arc types, with the exception of shortcuts, there exists a bi-directed version, where the gas can flow and the mentioned capabilities can be applied in both directions. While the mono-directed arcs only support positive flow w.r.t. their topological orientation, shortcuts are bi-directed by definition.
Auxiliary nodes do not have any special meaning and are used to decrease the number of necessary artificial arcs or to improve the general comprehensibility of the graph model. For the artificial arcs, we can often identify them with remotely controllable elements, which are contained in the original topology. But it is important to note that this does not hold in general. Combined arcs for example usually comprise at least one compressor station and one regulator.
Looking at the example in Fig. 3, we can identify the two anti-parallel regulators in the north-east of the original network topology with the bi-directed regulating arc in the station model. Additionally, we see two compressing arcs in the station model, which directly correspond to the two compressor stations depicted in the original topology. Here, the experts know that the compressor station in the east of the picture is used to compress gas coming from the north and leaving to the south, while the other, more central compressor station is used to compress gas coming from the south and leaving to the west. This explains their choice of endnodes and the direction of the artificial arcs.
Finally, the experts from OGE look at the possible operation modes and create the sets of possible flow directions and simple states. A flow direction consists of two subsets of fence nodes: Entries, where gas enters, and exits, where gas leaves. Additionally, they developed the simple states. Each simple state consists of a subset of flow directions that it supports and two subsets of artificial arcs: Arcs that have to be used and arcs that cannot be used. While an unusable arc can conceptually be seen as a closed valve, the former have to be used according to their corresponding models, which are described in Sect. 4.12. The goal for the design of the simple states is to summarize and approximate the technical capabilities at the original intersection area induced while keeping their overall number small.
In our example in Fig. 3, very common flow directions are situations, where gas enters from the north and/or east and leaves to the south and/or west. But as mentioned, it may also happen that gas enters from the south and leaves to the west. The creation of the corresponding simple states is then mainly based on experience and data regarding which of the operation modes of the corresponding intersection area have been used in the past.
4 Macroscopic transient gas flow model
In this section, we define our tri-level macroscopic transient gas flow model. We describe the entities of the underlying network and introduce variables and constraints representing their behaviour. Additionally, we explain the concepts describing their interplay and derive mathematical models for them.
Note that all constraints introduced here are part of the third level of our hierarchical MILP model, i.e., the level responsible for the technical control, and that the first and second levels do not feature any constraints. Furthermore, the first level does only control the inflow pressure slack variables, which are introduced in Sect. 4.3, and the second level is only in charge of the boundary value slack variables, which are introduced in Sect. 4.2. Both types of variables model the described non-technical measures, i.e., deviations from future pressure values at the sources and from supplies and demands, respectively. The assignment to the levels is due to the explained hierarchy and it is the goal of each level to minimize the extent of their usage, i.e., the sums of the absolute deviations, see Sect. 4.13 for the objectives and the overall model structure.
In the remainder of this paper, a gas network is modelled as a directed graph \(G = ({\mathcal {V}}, {\mathcal {A}})\), where \({\mathcal {V}} \) denotes the set of nodes and \({\mathcal {A}} \) the set of arcs.
4.1 Time steps and granularity
Additionally, we are given a set of time steps \({\mathcal {T}}_0 {:=} \{0, \dots , k \} \) together with a monotonically increasing function \(\tau : {\mathcal {T}}_0 \rightarrow {\mathbb {N}} \), called granularity. We assume that \(\tau (0) = 0\). In this context, \(\tau (t)\) represents the number of seconds that have passed until time step \(t \in {\mathcal {T}}_0 \) w.r.t. time step 0. For notational purposes let \({\mathcal {T}} {:=} {\mathcal {T}}_0 \setminus \{0\} \).
4.2 Boundary values
Furthermore, \({\mathcal {V}} ^{+} \subseteq {\mathcal {V}} \) and \({\mathcal {V}} ^{-} \subseteq {\mathcal {V}} \) denote the sources and sinks of the network, respectively, and we assume that \({\mathcal {V}} ^{+} \cap {\mathcal {V}} ^{-} = \emptyset \). While \({\mathcal {V}} ^\text {b} {:=} {\mathcal {V}} ^{+} \cup {\mathcal {V}} ^{-} \) is called the set of boundary nodes, \({\mathcal {V}} ^{0} {:=} {\mathcal {V}} \setminus {\mathcal {V}} ^\text {b} \) denotes the set of inner nodes.
For each boundary node \(v \in {\mathcal {V}} ^\text {b} \) and each time step \(t \in {\mathcal {T}} \) we are given so-called boundary values \(D_{v,t} \in {\mathbb {R}} \). They represent the future requirements in terms of supply (inflow), when \(v \in {\mathcal {V}} ^{+} \) is a source and we have \(D_{v,t} \in {\mathbb {R}}_{\ge 0} \), and demand (outflow), when \(v \in {\mathcal {V}} ^{-} \) is a sink and we have \(D_{v,t} \in {\mathbb {R}}_{\le 0} \). The boundary values may be adjusted dynamically to ensure the feasibility of the model. Thus, for each boundary node \(v \in {\mathcal {V}} ^\text {b} \) and \(t \in {\mathcal {T}} \) we introduce two continuous variables \(\sigma _{v,t}^{d+},\sigma _{v,t}^{d-} \in {\mathbb {R}}_{\ge 0} \). The actual boundary values, which are then considered in the model, are established through additional variables \(d_{v,t} \in {\mathbb {R}}_{\ge 0} \) for each source \(v \in {\mathcal {V}} ^{+} \) and \(d_{v,t} \in {\mathbb {R}}_{\le 0} \) for each sink \(v \in {\mathcal {V}} ^{-} \) and constraints
\(\sigma _{v,t}^{d+}\) and \(\sigma _{v,t}^{d-}\) are called boundary value slack variables and they are used to model the non-technical control measure of deviating from the given supplies and demands. Thus, they are controlled by the second level of our tri-level model, whose goal is to minimize their usage, i.e., their sum, see Sect. 4.13 for more details.
4.3 Pressures and pressure bounds
Additionally, for each node \(v \in {\mathcal {V}} \) we are given a non-negative pressure, which we denote by \(p _{v,0} \in {\mathbb {R}}_{\ge 0} \), representing the corresponding value in the initial state. Furthermore, we introduce pressure variables 
for each point in time \(t \in {\mathcal {T}} \). Here 
is a lower and \(\bar{p} _{v,t} \) is an upper bound on the pressure at node v and time t. These bounds are called technical pressure bounds.
For each boundary node \(v \in {\mathcal {V}} ^\text {b} \) and each point in time \(t \in {\mathcal {T}} \) we are additionally given so-called inflow pressure bounds
and
\(\bar{p} _{v,t} ^{\text {act}} \in {\mathbb {R}}_{\ge 0} \). These bounds are tighter than the technical pressure bounds and have to be respected if a boundary node has an nonzero boundary value. They represent the expected future pressures at the sources of the network. Nevertheless, in contrast to the hard technical pressure bounds they may be relaxed to ensure feasibility. Thus, we introduce two continuous variables
and
\(\sigma _{v,t}^{p-} \in [0, \bar{p} _{v,t} - \bar{p} _{v,t} ^{\text {act}} ]\) as well as constraints
\(\sigma _{v,t}^{p +}\) and \(\sigma _{v,t}^{p-}\) are called inflow pressure slack variables in the following and they are used to model the non-technical control measure of deviating from future source pressure values. Thus, they are controlled by the first level of our tri-level model, whose goal is to minimize their usage, i.e., their sum, see Sect. 4.13 for more details.
4.4 Mass flows
Next, we introduce variables representing the flow of gas on arcs in mass flow, which we are going to call simply flow in the following. Therefore, the arc set is partitioned into four sets \({\mathcal {A}} = {\mathcal {A}} ^{\text {va}} \, {\dot{\cup }} \, {\mathcal {A}} ^{\text {rg}} \, {\dot{\cup }} \, {\mathcal {A}} ^{\text {pi}} \, {\dot{\cup }} \, {\mathcal {A}} ^{\text {ar}} \), representing the different network element we consider. Here, \({\mathcal {A}} ^{\text {va}} \) denotes the set of valves, \({\mathcal {A}} ^{\text {rg}} \) the set of regulators (often synonymously called control valves in the literature), \({\mathcal {A}} ^{\text {pi}} \) the set of pipes, and \({\mathcal {A}} ^{\text {ar}} \) the set of so-called artificial arcs. The artificial arcs are further partitioned into mono-directed arcs \({\mathcal {A}} ^{\text {ar-mo}} \) and bi-directed arcs \({\mathcal {A}} ^{\text {ar-bi}} \), i.e., \({\mathcal {A}} ^{\text {ar}} = {\mathcal {A}} ^{\text {ar-mo}} \, {\dot{\cup }} \, {\mathcal {A}} ^{\text {ar-bi}} \), which is further discussed in Sect. 4.12. We allow parallel and anti-parallel arcs, but we do not allow loops.
For each mono-directed arc \(a \in {\mathcal {A}} ^{\text {rg}} \, \cup \, {\mathcal {A}} ^{\text {ar-mo}} \) and each time step \(t \in {\mathcal {T}} \) we introduce a variable \(q _{a,t} \in [0,\bar{q} _{a,t} ]\) representing the mass flow on the corresponding arc in forward direction. On the other hand, for valves and bi-directed artificial arcs we add two variables \(q _{a,t}^{\rightarrow },q _{a,t}^{\leftarrow } \in [0,\bar{q} _{a,t} ]\) representing mass flow in forward direction and backward direction on arc \(a \in {\mathcal {A}} ^{\text {va}} \,\cup \, {\mathcal {A}} ^{\text {ar-bi}} \), respectively. For pipes we distinguish in- and outflow to be able to account for changes in the amount of gas which is currently stored in the pipe. Therefore, for each pipe \(a = (\ell ,r) \in {\mathcal {A}} ^{\text {pi}} \) and each time step \(t \in {\mathcal {T}} \) we introduce two variables \(q _{\ell ,a,t},q _{r,a,t} \in [-\bar{q} _{a,t},\bar{q} _{a,t} ]\) representing the mass flow into a at \(\ell \) and out of a at r. Note that negative variable values represent mass flow out of a at \(\ell \) and into a at r, respectively. In all of the previous definitions, \(\bar{q} _{a,t} \) represents a practically reasonable flow bound.
Finally, for time step \(t = 0\) and each of the variables introduced above we are given an initial mass flow value, which is denoted with index 0.
4.5 Mass flow conservation
Next, for all nodes \(v \in {\mathcal {V}} \) we introduce mass flow conservation equations. For each inner node \(v \in {\mathcal {V}} ^{0} \) and each time step \(t \in {\mathcal {T}} \) the amount of flow entering v has to leave it and for a boundary node \(v \in {\mathcal {V}} ^\text {b} \) supply or demand must be satisfied. Hence, we have
For each inner node \(v \in {\mathcal {V}} ^{0} \) we introduce the same constraints except for the right side hand being 0.
4.6 Valves
Valves are network elements that can be used to link or unlink network parts by either creating a connection between the two corresponding endnodes or by disconnecting them. Thereby, a valve can be in one of two possible states. Either it is open, which implies that the pressure values at both ends are equal and mass flow is allowed in arbitrary direction (one can think of the endnodes being merged). Or it is closed, implying that there is no mass flow and the pressure values are independent or, as we synonymously call it, decoupled.
Thus, let \(a = (\ell ,r) \in {\mathcal {A}} ^{\text {va}} \) be a valve in \(G \). For each step in time \(t \in {\mathcal {T}} \) we introduce an additional binary variable \(z_{a,t} \in \{0,1\}\) indicating whether the valve is open or not. The behaviour described above can then be modelled using the following constraints
4.7 Regulators
Regulators can be seen as the continuous equivalent of valves. Besides being completely open or closed, regulators can also be partially open. Thereby, they generate friction due to which the gas pressure is decreased in the direction of flow. To model this behaviour, consider some \(a = (\ell ,r) \in {\mathcal {A}} ^{\text {rg}} \). For each \(t \in {\mathcal {T}} \) we add constraints
It is important to note that thereby the pressure at r can never be greater than the pressure at \(\ell \) in our model, i.e., we do not model so-called flap traps here. This mechanism closes the regulator if the pressure at r gets greater than the pressure at l and makes flow in the backward direction impossible. The reason for not including this mechanism in our model is that all regulators are considered to be connections to distribution parts of the network, i.e., parts only consisting of pipes, inner nodes and sinks, which are usually not operated at a higher pressure level than the upstream transportation network.
4.8 Pipes
One-dimensional gas flow in cylindric pipelines is usually described by the so-called Euler equations, a set of non-linear hyperbolic partial differential equations (Osiadacz 1996). In this paper we are going to assume isothermality, i.e., that the temperature remains constant. In that case the equations reduce to the Continuity Equation and the Momentum Equation. While the first equation ensures the conservation of mass, the second describes the interaction between the force acting on the gas particles and the rate of change in their momentum. Thus, for a pipe \(a = (\ell ,r) \in {\mathcal {A}} ^{\text {pi}} \) the isothermal Euler Equations can be stated as
Here, the x-variable represents the position in the pipe w.r.t. the distance from node \(\ell \). Furthermore, t denotes the time, and \(\rho \) and \(v \) the density and the velocity of the gas, respectively. Additionally, \(D_{a} \) denotes the diameter of the pipe and the gravitational acceleration is given by \(g \). Further, by \(\lambda _{a} \) we denote the friction factor of the pipe, which we derive from the formula of Nikuradse. This formula depends only on two characteristics of the pipe, namely its diameter and its integral roughness, see Fügenschuh et al. (2015) and Nikuradse (1950) for details. Finally, the slope of the pipe is \(s_{a} = \frac{h_r - h_\ell }{L_{a}} \in [-1,1]\), where \(h_\ell \) and \(h_r\) denote the altitude at \(\ell \) and r, respectively.
Next, we reformulate these equations w.r.t. the quantities we are interested in, i.e., mass flow \(q \), pressure \(p\), and velocity \(v\). First of all, mass flow is defined as
where \(A_{a} =D_{a} ^2\frac{\pi }{4}\) denotes the cross-sectional area of the pipe. Second, we apply the equation of state for real gases, which describes the relation between the gas pressure \(p\) and the density \(\rho \)
Here, \(z _{a}\) is the compressibility factor of the gas in the pipe and \(R _s\) denotes the specific gas constant. In the following, we assume both of them to be constants. For the compressibility factor this is a common assumption, see for example (Osiadacz 1996), and we define it as the average of the compressibility factors at both endnodes using their initial pressure values \(p _{\ell ,0} \) and \(p _{r,0} \) and the formula by Papay, see Saleh (2002). For the specific gas constant, this is a consequence of the fact that we assume the molar mass of the gas to be constant. Third, we drop the first and the third summand in the Momentum Equation, since their contribution under typical operating conditions in gas transport networks is negligible (Burlacu et al. 2019; Osiadacz 1996). Putting all this together, we can rewrite the equations and derive the so-called friction dominated model
Next, we discretize the equations using the implicit box scheme proposed in Domschke et al. (2011) and Kolb et al. (2010), respectively. Here, the length of the pipe \(L_{a} \) serves as spacial domain while we use the set of time steps \({\mathcal {T}}_0\) as time domain. Recall the definition of the mass flow variables \(q _{\ell ,a,t} \) and \(q _{r,a,t} \) from Sect. 4.5, where the first describes the flow into the pipe at node \(\ell \) and the second the flow out of the pipe at node r. The discretized equations for two adjacent time steps \(t-1\) and t, where \(t \in {\mathcal {T}} \), can then be written as
Finally, we are going to use the linear model for the Momentum Equation, which was proposed by Hennings (2018). To derive it, we fix the absolute velocities in the friction-based pressure difference term of the Momentum Equation, i.e., in the third summand, to the absolute gas velocities in the initial time step. Thereby, we derive the following equations, which we use to model the transient gas flow in pipelines
The obvious issue arising here is that if the velocity of the mass flow in- or decreases significantly over time, we might under- or overestimate the friction loss, respectively. An analysis of historic real-world data in Hennings (2018) shows, that the resulting deviation can be significant, but is small in most of the cases when considering a time horizon of 12 h. However, in Sect. 6 we introduce and discuss an iterative velocity adjustment procedure inspired by sequential linear programming, which we apply in order to determine solutions that are also feasible for the non-linear Momentum Equation (M).
4.9 Network stations
The main idea behind the network station model and the process of how it is derived are discussed in Sect. 3. Formally, within \(G \) there exist \(m \in {\mathbb {N}} \) subgraphs \(G_{i} =({\mathcal {V}} _{i},{\mathcal {A}} ^{\text {ar}}_{i}) \) called network stations, which consist of inner nodes and artificial arcs only, i.e., \({\mathcal {V}} _{i} \subseteq {\mathcal {V}} ^{0} \) and \({\mathcal {A}} ^{\text {ar}}_{i} \subseteq {\mathcal {A}} ^{\text {ar}} \) for all \(i \in \{1, \dots , m \}\). Each artificial arc is contained in exactly one station and each inner node is contained in at most one station, i.e., \({\mathcal {A}} ^{\text {ar}}_{i} \cap {\mathcal {A}} ^{\text {ar}}_{j} = \emptyset \) and \({\mathcal {V}} _{i} \cap {\mathcal {V}} _{j} = \emptyset \) holds for \(i,j \in \{1,\ldots ,m \}\) with \(i \ne j\) and we have \({\mathcal {A}} ^{\text {ar}} = \bigcup \nolimits _{i = 1}^{m} {\mathcal {A}} ^{\text {ar}}_{i} \).
The node set \({\mathcal {V}} _{i} \) can be further partitioned into so-called fence nodes \({\mathcal {V}} ^{\text {fn}} _{i} \) and auxiliary nodes \({\mathcal {V}} ^{\text {ar}} _{i} \), i.e., \({\mathcal {V}} _{i} {:=} {\mathcal {V}} ^{\text {fn}} _{i} \, {\dot{\cup }} \, {\mathcal {V}} ^{\text {ar}} _{i} \). A node \(v \in {\mathcal {V}} _{i} \) is a fence node if it is connected to at least one arc outside the gas network station, i.e., if \(\delta (v) \cap ({\mathcal {A}} ^{\text {pi}} \, \cup \, {\mathcal {A}} ^{\text {rg}} \, \cup \, {\mathcal {A}} ^{\text {va}}) \ne \emptyset \), where \(\delta (v)\) denotes the set of arcs incident to v. Otherwise, if \(\delta (v) \subseteq {\mathcal {A}} ^{\text {ar}}_{i} \), it is an auxiliary node.
Additionally, \({\mathcal {F}} _{i} \subseteq {\mathcal {P}}({\mathcal {V}} ^{\text {fn}} _{i}) \times {\mathcal {P}}({\mathcal {V}} ^{\text {fn}} _{i})\) denotes the set of so-called flow directions of gas network station \(G_{i} \), where \({\mathcal {P}}\) is the powerset operator. A flow direction \(f = (f ^{+},f ^{-}) \in {\mathcal {F}} _{i} \) consists of its entry fence nodes \(f ^{+} \subseteq {\mathcal {V}} ^{\text {fn}} _{i} \) and its exit fence nodes \(f ^{-} \subseteq {\mathcal {V}} ^{\text {fn}} _{i} \) and it holds that \(f ^{+} \cap f ^{-} = \emptyset \).
Furthermore, the set \({\mathcal {S}}_{i} \subseteq {\mathcal {P}}( {\mathcal {F}} _{i}) \times {\mathcal {P}}({\mathcal {A}} ^{\text {ar}}_{i}) \times {\mathcal {P}}({\mathcal {A}} ^{\text {ar}}_{i})\) containing the so-called simple states is given for each gas network station \(G_{i} \). A single simple state \(s = (s ^{f}, s ^{\text {on}}, s ^{\text {off}}) \in {\mathcal {S}}_{i} \) is composed of the set of flow directions \(s ^{f} \) it supports as well as the set of its active \(s ^{\text {on}} \) and its inactive artificial arcs \(s ^{\text {off}} \).
4.10 Controlling network stations
In each time step \(t \in {\mathcal {T}}_0 {:=} \{0, \dots , k \} \) three types of control decisions have to be made for a gas network station \(G_{i} \). These decisions impact each other and can be put into a hierarchical order. Here, we describe an order in a top to bottom fashion and introduce the variables and constraints modelling the decisions and their interplay.
First of all, exactly one flow direction \(f \in {\mathcal {F}} _{i} \) has to be chosen for each \(G_{i} \). Given this flow direction, one must additionally choose exactly one simple state \(s \in {\mathcal {S}}_{i} \) which supports this flow direction, i.e., \(f \in s ^{f} \) has to hold. Given a decision on the simple state, all arcs in \(s ^{\text {on}} \) must be active, while the inactive arcs \(s ^{\text {off}} \) cannot be used. For all remaining artificial arcs \(a \in {\mathcal {A}} ^{\text {ar}}_{i} \setminus (s ^{\text {on}} \cup s ^{\text {off}})\), which we call optional arcs, we can independently choose whether they are active or not.
Thus, for each time step \(t \in {\mathcal {T}}_0 \) we introduce binary variables \(x_{f,t} \in \{0,1\}\) for each flow direction \(f \in {\mathcal {F}} _{i} \), \(x_{s,t} \in \{0,1\}\) for each simple state \(s \in {\mathcal {S}}_{i} \), as well as \(x_{a,t} \in \{0,1\}\) for each artificial link \(a \in {\mathcal {A}} ^{\text {ar}}_{i} \) all indicating whether the corresponding entity is active at that point in time or not. Furthermore, for each network station \(G_{i} \) we add the following constraints
While constraints (12) and (13) ensure that exactly one flow direction and one simple state are chosen for each time step \(t \in {\mathcal {T}}_0 \), (14) guarantees that the chosen simple state supports the chosen flow direction. Additionally, constraints (15) and (16) make sure that the artificial arcs corresponding to the simple state are active or not, respectively. No condition is imposed on the optional arcs.
Next, in order to penalize changes over time w.r.t. simple states or artificial links in the objective function we introduce additional binary variables. For each station \(G_{i} \) and each time step \(t \in {\mathcal {T}} \) we have \(\delta _{s,t} \in \{0,1\}\) for each \(s \in {\mathcal {S}}_{i} \) and \(\delta _{a,t} ^{\text {on}},\delta _{a,t} ^{\text {off}} \in \{0,1\}\) for each \(a \in {\mathcal {A}} ^{\text {ar}}_{i} \). Furthermore, we add constraints
While \(\delta _{s,t} \) and \( \delta _{a,t} ^{\text {on}} \) indicate whether or not a simple state or artificial link has been switched on in time step t, \(\delta _{a,t} ^{\text {off}} \) additionally indicates whether or not an artificial link has been switched off. For the simple states we do not need such a variable, since we know that exactly one of them is active in each time step, but in the case of optional artificial arcs this does not apply. All variables \(\delta _{s,t} \) are associated with an individual cost parameter \(w ^s \in {\mathbb {R}}_{\ge 0} \), while variables \(\delta _{a,t} ^{\text {on}} \) as well as \(\delta _{a,t} ^{\text {off}} \) are assigned a common cost parameter \(w ^a \in {\mathbb {R}}_{\ge 0} \).
Note that we refrain from penalizing changes w.r.t. the flow directions here. For example, if a network station is in a so-called bypass state, e.g., only shortcuts are active, the gas may “slosh” back and forth between the fence nodes. Thus, flow directions changes do not suggest an unstable behaviour and are actually very common.
4.11 Flow direction related constraints
Activating a flow direction imposes certain conditions on the mass flow and pressure values w.r.t. a gas network station \(G_{i} \). Most importantly, for a flow direction \(f = (f ^{+},f ^{-}) \in {\mathcal {F}} _{i} \) no outflow is allowed at its entry and no inflow at its exit fence nodes. It is however allowed that there is no flow at all, which is the condition that must hold for all other fence nodes \(v \in {\mathcal {V}} ^{\text {fn}} _{i} \setminus (f ^{+} \cup f ^{-})\). Furthermore, for some of the fence nodes there exist additional pressure bounds if a flow direction is chosen in which they serve as exits. And finally, there exist conditions on the sums of absolute amounts of flow of subsets of fence nodes, which have to be satisfied in order to activate certain flow directions.
4.11.1 In- and outflow constraints
First of all, for each fence node \(v \in {\mathcal {V}} _{i} \) and each point in time \(t \in {\mathcal {T}}_0 \) we introduce two continuous variables \(q _{v,t} ^\text {in},q _{v,t} ^\text {out} \in {\mathbb {R}}_{\ge 0} \) that, together with the following constraint, account for the total in- or outflow from outside the station, respectively
Note that one could alternatively sum up the mass flow values of the incident pipes, regulators, and valves on the left hand side and switch the signs of the variables on the right hand side of the equation. This is because flow conservation holds at the fence nodes, since \({\mathcal {V}} ^{\text {fn}} _{i} \subseteq {\mathcal {V}} ^{0} \). Next, for each flow direction \(f = (f ^{+},f ^{-}) \in {\mathcal {F}} _{i} \) we introduce the following constraints:
Here, \({\overline{q}}^{\text {in}}_{v,t}\) and \({\overline{q}}^{\text {out}}_{v,t}\) are upper and lower bounds on the maximum possible in- and outflow, respectively, which can be derived from constraints (19) together with the above mentioned alternative constraint. If a flow direction is active, \(q _{v,t} ^\text {in} \) can be nonzero for the entry and \(q _{v,t} ^\text {out} \) for the exit fence groups only.
4.11.2 Exit pressure bounds
Furthermore, for some fence nodes \(v \in {\mathcal {V}} ^{\text {fn}} _{i} \) there exists an additional upper pressure bound \(\bar{p} _{v}^{\text {exit}} \), which is tighter than its technical upper bound and has to be respected if a flow direction \(f = (f ^{+},f ^{-}) \in {\mathcal {F}} _{i} \) is active, for which v is an exit fence node, i.e., for which \(v \in f ^{-} \). This can be modelled via the following constraints
4.11.3 Flow direction conditions
Finally, for each network station, there exists a set of so-called flow direction conditions \({\mathcal {W}} _{i} \subseteq {\mathcal {F}} _{i} \times {\mathcal {P}} ({\mathcal {V}} ^{\text {fn}} _{i}) \times {\mathcal {P}} ({\mathcal {V}} ^{\text {fn}} _{i})\), demanding that the sum of the absolute in- and outflows of the first set of fence nodes is less than or equal than the sum of the in- and outflows of the second set in order to activate the corresponding flow direction. Hence, for each \(w = (f,{\mathcal {V}} _{w_1},{\mathcal {V}} _{w_2}) \in {\mathcal {W}} _{i} \) we introduce
where \(M_{w,t} {:=} \sum \nolimits _{v \in {\mathcal {V}} _{w_1} \cap f ^{+}} {\overline{q}}^{\text {in}}_{v,t} + \sum \nolimits _{v \in {\mathcal {V}} _{w_1} \cap f ^{-}} {\overline{q}}^{\text {out}}_{v,t}\). They are introduced because certain simple states with the ability to compress need these conditions in order to work.
4.12 Artificial arcs
The set of artificial arcs can be further partitioned into four disjoint subsets \({\mathcal {A}} ^{\text {ar}} = {\mathcal {A}} ^{\text {ar-sc}} {\dot{\cup }} {\mathcal {A}} ^{\text {ar-rg}} {\dot{\cup }} {\mathcal {A}} ^{\text {ar-co}} {\dot{\cup }} {\mathcal {A}} ^{\text {ar-cb}} \). Here, \({\mathcal {A}} ^{\text {ar-sc}} \) denotes the set of so-called shortcuts, \({\mathcal {A}} ^{\text {ar-rg}} \) the set of so-called regulating arcs, \({\mathcal {A}} ^{\text {ar-co}} \) the set of so-called compressor arcs, and \({\mathcal {A}} ^{\text {ar-cb}} \) the set of so-called combined arcs. Further, we denote the set of pressure increasing arcs by \({\mathcal {A}} ^{\text {ar-pr}} = {\mathcal {A}} ^{\text {ar-co}} \cup {\mathcal {A}} ^{\text {ar-cb}} \). The sets \({\mathcal {A}} ^{\text {ar-sc}} _{i} \subseteq {\mathcal {A}} ^{\text {ar-sc}} \), \({\mathcal {A}} ^{\text {ar-rg}} _{i} \subseteq {\mathcal {A}} ^{\text {ar-rg}} \), \({\mathcal {A}} ^{\text {ar-co}} _{i} \subseteq {\mathcal {A}} ^{\text {ar-co}} \), \({\mathcal {A}} ^{\text {ar-cb}} _{i} \subseteq {\mathcal {A}} ^{\text {ar-cb}} \), and \({\mathcal {A}} ^{\text {ar-pr}} _{i} \subseteq {\mathcal {A}} ^{\text {ar-pr}} \) describe the corresponding entities contained in gas network station \(G_{i} \).
In this section we explain how the artificial arcs and their different capabilities when controlling the gas flow through the station are modelled. But first, we shortly explain the difference between bi-directed and mono-directed arcs.
4.12.1 Bi-Directed arcs
In contrast to the mono-directed arcs, mass flow and pressure modifications according to the corresponding artificial arc are possible into both directions on bi-directed arcs. Further, their capabilities, for example the compression of gas, are also applicable in both directions. Thus, in our model for these arcs we first decide into which direction the mass flow is going at each point in time.
Therefore, for each bi-directed arc \(a \in {\mathcal {A}} ^{\text {ar-bi}} \) and each time step \(t \in {\mathcal {T}}_0 \) we introduce two binary variables \(x_{a,t} ^{\rightarrow },x_{a,t} ^{\leftarrow } \in \{0,1\}\) encoding the direction of the flow in case the arc is active and add constraints
to the model. Given this decision, bi-directed arcs are modelled analogously to their mono-directed counterparts using the corresponding variable.
4.12.2 Shortcuts
All shortcuts are bi-directed arcs and mass flow is possible into both directions. They can conceptually be seen as the equivalent of valves (see Sect. 4.6) inside a station and are used to connect network parts if the corresponding pressure levels are equal. Thus, for each shortcut \(a = (\ell ,r) \in {\mathcal {A}} ^{\text {ar-sc}} \) we add constraints
If a shortcut is active at time \(t \in {\mathcal {T}} \), i.e., if \(x_{a,t} = 1\), the pressures at \(\ell \) and r have to be equal and mass flow can go into an arbitrary direction with an arbitrary value, i.e., there may be be forward flow \(q _{a,t}^{\rightarrow } \in [0,\bar{q} _{a,t} ]\) or backward flow \(q _{a,t}^{\leftarrow } \in [0,\bar{q} _{a,t} ]\) depending on the decision made in constraint (24). If the shortcut is not active, the pressure values are decoupled and there is no flow.
4.12.3 Regulating arcs
Regulating arcs can conceptually be seen as the equivalent of regulators (see Sect. 4.7) inside a gas network station. They are used to decrease the gas pressure in the direction of the mass flow, which is needed if, for example, gas enters a distribution network, which is technically not suited for the high pressure of the transportation network. Hence, for regulating arcs \(a = (\ell ,r) \in {\mathcal {A}} ^{\text {ar-rg}} \) we introduce the following constraints
If a regulating arc is active at some point in time \(t \in {\mathcal {T}} \), i.e., if \(x_{a,t} = 1\), the pressure at \(\ell \) has to be greater or equal than the pressure at r. Otherwise, the pressure values are decoupled and there is no mass flow. For bi-directed regulating arcs \(a = (\ell ,r) \in {\mathcal {A}} ^{\text {ar-rg}} \cap {\mathcal {A}} ^{\text {ar-bi}} \), we derive an analogous set of constraints using \(x_{a,t} ^{\rightarrow } \) and \(x_{a,t} ^{\leftarrow } \) instead of \(x_{a,t} \).
4.12.4 Pressure increasing arcs
The pressure increasing arcs \({\mathcal {A}} ^{\text {ar-pr}} \), i.e., the compressor arcs \({\mathcal {A}} ^{\text {ar-co}} \) and the combined arcs \({\mathcal {A}} ^{\text {ar-cb}} \), are key elements when it comes to control a macroscopic gas network. They are able to compress gas and thereby increase its pressure, which makes up for pressure loss due to friction in the pipes or height differences that have to be overcome.
In our model, one can conceptually think of one (big) compressor unit being installed at each arc \(a \in {\mathcal {A}} ^{\text {ar-pr}} _{i} \) of each gas network station \(G_{i} \). The maximum power it has available for compression \({\tilde{\pi }}_{a,t} \in {\mathbb {R}}_{\ge 0} \), the maximum amount of mass flow that can pass through it \({\tilde{q}}_{a,t} \in {\mathbb {R}}_{\ge 0} \), and its maximum compression ratio \({\tilde{r}}_{a,t} \in [1,\infty )\) are dynamically determined in each time step through an assignment of approximations of real-world compressor units, simply called machines in the following, and a linear combination of their corresponding values.
Thus, for each station \(G_{i} \) we are given a set of machines \({\mathcal {M}}_{i} \) and each machine \(m \in {\mathcal {M}}_{i} \) possesses an associated power value \(P _{m,t} \in {\mathbb {R}}_{\ge 0} \), a maximum mass flow \(Q _{m,t} \in {\mathbb {R}}_{\ge 0} \), and a maximum compression ratio \(R_{m,t} > 1\) for each time step \(t \in {\mathcal {T}} \). Further, for each pressure increasing arc \(a \in {\mathcal {A}} ^{\text {ar-pr}} _{i} \) there exists a subset of machines \({\mathcal {M}}_{i} ^{a} \subseteq {\mathcal {M}}_{i} \) that can potentially be assigned to it and a maximum number of assignable machines \(M_a^{\text {max}} \). Since each machine can be assigned to at most one compressing link in each time step \(t \in {\mathcal {T}} \), we introduce binary variables \(y_{m,a,t} \in \{0,1\}\) indicating whether machine \(m \in {\mathcal {M}}_{i} \) is assigned to arc \(a \in {\mathcal {A}} ^{\text {ar-pr}} _{i} \) or not, and add constraints
In the real world, compressor units are operated alone, in parallel, sequentially or in a parallel-sequential setting. This is achieved by the opening and closing of valves in the surrounding piping. By setting them up in parallel, a larger amount of mass flow can be compressed, while in serial a higher compression ratio can be achieved. In our model, we refrain from choosing a setup for the machines and overestimate the capabilities of pressure increasing arcs in the sense that we assume that the maximum amount of flow (parallel setting) and the highest compression ratio (sequential setting) are available at the same time. Thus, we add the following constraints
The first constraint (33) determines the power available on arc \(a \in {\mathcal {A}} ^{\text {ar-pr}} \) by adding up the maximum power of the assigned machines. Analogously, the second constraint (34) determines the maximum amount of mass flow that can be compressed. On the other hand, the third constraint (35) is a (conservative) approximation of the maximum compression ratio, which is used in order to avoid non-linear constraints.
Finally, the connection between pressure difference, the amount of mass flow passing through a compressor machine, and the power necessary to realize it is given by the non-linear power equation for compressor machines
where \(\pi _{a,t} \in {\mathbb {R}}_{\ge 0} \) is the variable representing the necessary power when a mass flow of \(q _{a,t} \) with initial pressure \(p _{\ell ,t} \) shall be compressed up to \(p _{r,t} \). Here, \(\eta _\text {ad} \) is the adiabatic efficiency of the compression, which we assume to be constant for all existing compressor machines, and \(\kappa = 1.296\) [14] is the isentropic exponent.
To avoid introducing this non-linear constraint we determine a linear approximation as follows: For each artificial compressing link \(a \in {\mathcal {A}} ^{\text {ar-pr}} \) and each \(t \in {\mathcal {T}} \), we sample N points 
, where \(\pi ^\text {max}_{a,t}\) is the maximum possible power for a at t derived from (31) and (32), such that \(p _{\ell ,t} \le p _{r,t} \) and determine the corresponding mass flow \(q _{a,t} \) using the original power equation. To the resulting set of 4-tuples we apply an ordinary least-squares method and determine coefficients \((\alpha _0,\alpha _1,\alpha _2,\alpha _3)\) for a linear approximation, which gives rise to constraints
where we assume that \(\alpha _1 \in {\mathbb {R}}_{\le 0} \) and \(\alpha _2 \in {\mathbb {R}}_{\ge 0} \) (otherwise we use the corresponding other bound for the coefficients of \(x_{a,t} \) on the right hand sides). If the pressure increasing arc is active, it has to respect this linear approximation. Otherwise, there is no flow and the pressures at both ends are decoupled. Here, typical values of the coefficients are in the ranges \(\alpha _0 \in [500,25000]\), \(\alpha _1 \in [-70,-10]\), \(\alpha _2 \in [3,40]\), and \(\alpha _3 \in [400,1200]\). However, the concrete values depend on the number of assignable machines, their properties, as well as different mass flow and pressure bounds. To give a concrete example, for an active compressing arc (\(x_{a,t} \) = 1) which features up to four assignable machines, we derive the following relationship
when constraints (36) and (37) are merged.
Finally, we add the following set of constraints
The first two constraints (38) and (39) ensure that the mass flow and power used for compression do not violate the upper bounds given by the machine assignments. Finally, the outgoing pressure is bounded by the product of the initial ingoing pressure at \(t=0\) and the current maximum compression ratio (40) if the corresponding arc is active. Using the pressure variables instead would again result in non-linear constraints.
4.12.5 Compressor arcs
Besides constraints (31)–(40), for each compressor arc \(a = (\ell ,r) \in {\mathcal {A}} ^{\text {ar-co}} \) and each time step \(t \in {\mathcal {T}} \) we add constraints
If the arc is active at some point in time \(t \in {\mathcal {T}} \), i.e., \(x_{a,t} = 1\), the pressure at \(\ell \) has to be smaller than or equal to the pressure at r. Further, we bound \(p _{r,t} \) by \(r^{\text {max}}_{a,t} p _{\ell ,t} \) where \(r^{\text {max}}_{a,t} \) is the maximum possible compression ratio of a at time t, which can be derived from constraints (32) and (40). If it is not active, the pressure values are decoupled and there is no mass flow due to constraints (39) and (34).
Further, there may be an additional upper bound on the pressure at node r, which has to be respected if the arc is active. Let \(\bar{p} _{r,t} ^{\text {out}} \) denote this upper bound. We can model this requirement by:
If such a bound is given, we shrink the sample space described in the previous section, accordingly.
4.12.6 Combined arcs
A combined arc \(a=(\ell ,r) \in {\mathcal {A}} ^{\text {ar-cb}} _{i} \) can be used as a regulating or a compressor arc. Hence, we first of all introduce two binary decision variables encoding in which mode it is activated.
All constraints (31)–(40), where \(x_{a,t} \) is replaced by \(x_{a,t} ^{\text {cp}} \) in (32), (36), (37), and (40), are added for each combined arcs except for (39), which is replaced by
since \({\tilde{q}}_{a,t} = 0\) holds if \(x_{a,t} ^{\text {rg}} = 1\). To capture the behaviour as regulating arc, we add constraints
analogously to (9), while for the pressure increasing arc we additionally have
analogously to (41) and (42). Further, as for the compressor arcs, there may be an additional upper bound on the pressure at node r, if compression is used. Thus,
is added, similar to (43), and we shrink the sample space from the previous section, accordingly.
4.13 Objectives and complete model
As mentioned before, we solve a hierarchical MILP formulation consisting of three levels, i.e., a tri-level mixed-integer linear program. This is motivated by the following rationale. In the real-world dispatchers first of all try to control the network by using technical measures only, i.e., by using and changing the setting of the remotely controllable elements in order to satisfy the supplies and demands desired by the customers. If this does not seem to work, they have some non-technical measures at hand. The most common and standardized ones are the change of supplies and demands by either buying or selling gas, i.e., using so-called balancing energy, or by using contractual options like the interruption of customers. If changing the supplies and demands alone does not work out, the last option is to ask other transport system operators for pressure changes of the future supply at some source nodes. In practice, this is done by phone calls and can therefore be seen as the last possible non-standardized option. Thus, we finally state the complete tri-level MILP formulation as
In our model the first level controls the slack variables for the inflow pressure bounds, while the second level controls the slack variables for the boundary values. The goal of both levels is to minimize the sum of the corresponding slack variables, i.e., the sum of absolute deviations. The third level, which can be seen as the level for the technical operation of the network while the other two only ensure feasibility, controls the remaining variables and minimizes the total cost, which is given as the weighted sum of simple state and auxiliary link changes.
5 An algorithmic approach for solving the tri-level MILP model
To state our algorithmic approach for solving the tri-level MILP model, which we call AATM in the following, we define three closely connected (single-level) MILP formulations: First, the MILP consisting of objective function (52) and constraint set (1)–(49) with all slack variables being fixed to zero, which we denote by \(L_3\) in the following. Second, objective function (51) together with (1)–(49) as well as all inflow pressure slack variables being fixed to zero we denote as MILP \(L_2\). And third, (50) combined with (1)–(49) describes MILP formulation \(L_1\), which is often synonymously called high point relaxation in the context of hierarchical optimization. The AATM is stated in Algorithm 1.
Within the AATM we may use any exact algorithm, which can prove whether the MILPs defined above are infeasible or determine an optimal solution for it. Note that if there exists a feasible solution for any of them, there exists an optimal solution since the objective functions are all bounded from below by 0. Typically, linear-programming (LP) based branch-and-bound algorithms are applied, such as the one provided by Gurobi Optimization (2020), which we used for our computational experiments in Sect. 8. For more information, we refer to Achterberg (2007).
If there exists a feasible solution with no slacks, i.e., for \(\text {L}_3\), an optimal solution for it is returned as an optimal solution for L in line (23). Otherwise, if there exists a feasible solution for \(\text {L}_2\), we subsequently solve \(\hat{\text {L}}_3\), i.e., \(\text {L}_3\) with all slack variables fixed to an optimal solution of \(\text {L}_2\). Doing this, we again determine an optimal solution for the tri-level MILP, see lines (18)–(20). Finally, if \(\text {L}_2\) does not admit a feasible solution, we consider the high point relaxation \(\text {L}_1\). If it is infeasible, the whole problem L is infeasible (7). Otherwise, we subsequently solve \(\tilde{\text {L}}_2\) and \(\tilde{\text {L}}_3\) to determine an optimal solution for L, see lines (9)–(15).
5.1 Heuristics for single-level MILP formulations
Next, we introduce two heuristics that can be applied to the single-level MILP formulations from the previous section, i.e., to \(\text {L}_1, \text {L}_2, \tilde{\text {L}}_2, \text {L}_3, \tilde{\text {L}}_3,\) or \(\hat{\text {L}}_3\). Having branch-and-bound based algorithms in particular in mind, our main motivation is to produce feasible solutions, use them within the chosen exact algorithm to solve MILPs in Algorithm 1, and thereby try to accelerate the solving process.
First, we introduce a rolling horizon approach. We start by solving the model for the variables and constraints corresponding to \(t = 0\) only. Next, in each of the following n iterations we solve the model where the next time step is added and the binary decisions of all but the newly added time step are fixed to the solution values from the previous iteration.
The second heuristic initially solves a specially designed Min-Cost-Flow (MCF) problem defined on \(G = ({\mathcal {V}}, {\mathcal {A}})\) for each time step \(t \in {\mathcal {T}} \). Analyzing the in- and outflows at the fence nodes of each station in the optimal solutions, we reduce the number of possible flow directions by fixing binary variables of those flow directions for this time step to zero, which are not consistent with the MCF solution.
5.1.1 Rolling horizon heuristic
The first idea to determine a feasible solution for any of the MILP models above is to use a rolling horizon approach, i.e., to iteratively consider the model with an additional time step, solve it, and fix the binary decision variables in the next iteration to the solution values of the corresponding binary variables of an optimal solution of the current iteration. This procedure is stated in Algorithm 2.
Let \(\text {MILP}_{k}^{j}\) denote the model on the first k time steps with all binary variables corresponding to some \(i \in {\mathcal {T}}_0 \) with \(i < j\) fixed. If a MILP formulation \(\text {MILP}_{k}^{k}\) turns out to be infeasible, we stop the heuristic and return UNSUCCESSFUL. Otherwise, the heuristic terminates with a feasible solution. Rolling horizon style approaches have been widely used to find feasible solutions for time-dependent optimization problems, for example for disruption management in the railway industry (Nielsen et al. 2012) or scheduling problems (Addis et al. 2016; Samà et al. 2013).
5.1.2 A min-cost-flow based heuristic
The idea behind the MCF heuristic is to decrease the number of binary variables regarding the flow directions in the MILP formulation by fixing a subset of them for each time step to zero. Due to the hierarchical structure of the decisions in a station, i.e., the constraints modelling the interplay between the binary variables corresponding to flow directions, simple states, and artificial arcs, further fixations may follow. In order to exclude certain flow directions, we solve a specially designed MCF problem on the underlying graph for each time step and analyze the in- and outflows at the fence nodes. Besides the existence of fast algorithms for solving it, MCF is often used in practice by TSOs to approximate the flow of natural gas, see for example (Hoppmann and Schwarz 2018; Steringa et al. 2015), where it is used to determine worst case transport scenarios.
For each arc \(a \in {\mathcal {A}} ^{\text {rg}} \cup {\mathcal {A}} ^{\text {ar-mo}} \) we introduce a non-negative flow variable \(f _{a} \in {\mathbb {R}}_{\ge 0} \) and for each \(a \in {\mathcal {A}} ^{\text {va}} \cup {\mathcal {A}} ^{\text {pi}} \cup {\mathcal {A}} ^{\text {ar-bi}} \) we introduce two non-negative flow variables \(f _{a}^{\rightarrow },f _{a}^{\leftarrow } \in {\mathbb {R}}_{\ge 0} \) describing the forward and the backward flow, respectively. Further, consider some \(t \in {\mathcal {T}} \) and w.l.o.g. we have \(\sum \nolimits _{v \in {\mathcal {V}} ^{+}} D_{v,t} \ge \sum \nolimits _{v \in {\mathcal {V}} ^{-}} |D_{v,t} | > 0\) and let \(\chi _t {:=} \frac{\sum \nolimits _{v \in {\mathcal {V}} ^{+}} D_{v,t} }{\sum \nolimits _{v \in {\mathcal {V}} ^{-}} |D_{v,t} |}\). Additionally, recall that for each pipe \(a \in {\mathcal {A}} \) we are given its length \(L_{a} \in {\mathbb {R}}_{\ge 0} \). Finally, the MCF we solve for each time step t can be formulated as the following linear program (LP)
Note that if there is no supply or demand in a time step, we define the right hand sides of all constraints to be 0. Now, given an optimal solution, for each gas network station \(G_{i} \) and each of its fence nodes \(v \in {\mathcal {V}} ^{\text {fn}} _{i} \) we check whether there is in- or outflow w.r.t. to \(G_{i} \) at v. If for
we have \(f_v \in {\mathbb {R}}_{\ge 0} \setminus [0,\gamma )\) we call v an MCF entry fence node and if \(f_v \in {\mathbb {R}}_{\le 0} \setminus (\gamma ,0]\) we call it an MCF exit fence node. The idea to use some threshold \(\gamma \in {\mathbb {R}}_{\ge 0} \) in this definition is that small in- and outflows could be realized by using the gas stored in adjacent pipelines, i.e., linepack.
Next, we call a flow direction \(f = (f ^{+},f ^{-})\) MCF-valid for station \(G_{i} \) and time step t, if for each MCF entry fence node v we have \(v \in f ^{+} \) and for each MCF exit fence node v we have \(v \in f ^{-} \). As the final step of the heuristic, we solve the MILP model as described in Sect. 4, where we fix for each station \(G_{i} \) and each time step \(t \in {\mathcal {T}} \) all binary variables corresponding to flow directions which are not MCF-valid to zero, if there exists at least one MCF-valid flow direction. Otherwise, we do not apply any fixations.
5.2 Solution smoothing
Due to the nature of LP-based branch-and-bound algorithms, such as the one provided by Gurobi (2020) which we use for our computational experiments in Sect. 8, a typical issue that arises is the non-smoothness of solutions. This can for example be observed on compressing arcs, where in many cases massive amounts of gas are compressed in a single time step, while there is basically no compression at all in all other time steps. Instead, we would like to have a steady mass flow over the whole time horizon, if that is possible. Similarly, the same behaviour can be observed for the outgoing pressures. Of course, these solutions with big differences in the corresponding variable values for consecutive time steps may be feasible w.r.t. the model, but such a behaviour is not desirable in practice. In the following, we consider a solution to be smooth if the changes in the inflow and pressure values at the fence nodes of the network stations are as small as possible for two consecutive time steps. Thus, we introduce the following LP-based smoothing routine.
Given a solution S for the hierarchical MILP formulation, consider the following single-level LP: First, fix all binary variables to their corresponding solution values. If all inflow pressure slack variables are zero, fix them as well. If furthermore all supply and demand slacks are zero, we fix them, too. We denote this induced linear program by \(\text {LP}_{\text {base}}(S)\) in the following.
Next, for each network station \(G_{i} \), each fence node \(v \in {\mathcal {V}} ^{\text {fn}} _{i} \), and each time step \(t \in {\mathcal {T}} \) we add four continuous variables \(\delta _{v,t}^{p +}, \delta _{v,t}^{p-}, \delta _{v,t}^{q +}, \delta _{v,t}^{q-} \in {\mathbb {R}}_{\ge 0} \) as well as two continuous variables \({\overline{\delta }}_{v}^{p}, {\overline{\delta }}_{v}^{q} \in {\mathbb {R}}_{\ge 0} \) together with the following constraints:
Constraint (53) measures the difference between the pressure values at v from time step \(t-1\) to time step t using variables \(\delta _{v,t}^{p +}\) and \(\delta _{v,t}^{p-}\). Analogously, the difference between in- and outflow is measured by constraint (55) using variables \(\delta _{v,t}^{q +}\) and \(\delta _{v,t}^{q-}\). The maximum difference between any two consecutive time steps w.r.t. pressure and in- and outflow at v is determined by constraints (54) and (56) using variables \({\overline{\delta }}_{v}^{p}\) and \({\overline{\delta }}_{v}^{q}\), respectively.
Furthermore, for each network station \(G_{i} \) and each time step \(t \in {\mathcal {T}} \) we add two continuous variables \({\overline{\delta }}_{i,t}^{p}, {\overline{\delta }}_{i,t}^{q} \in {\mathbb {R}}_{\ge 0} \) and introduce constraints
Here, we determine the maximum difference w.r.t. to pressure and in- and outflow between time step \(t-1\) and time step t at the fence nodes of station \(G_{i} \) using \({\overline{\delta }}_{i,t}^{p}\) and \({\overline{\delta }}_{i,t}^{q}\), respectively.
While variables \({\overline{\delta }}_{v}^{p}\) and \({\overline{\delta }}_{v}^{q}\) are associated with objective coefficients \(w _{v}^{\text {sm-q}} \) and \(w _{v}^{\text {sm-p}} \), \({\overline{\delta }}_{i,t}^{p}\) and \({\overline{\delta }}_{i,t}^{q}\) are associated with objective coefficients \(w _{i}^{\text {sm-q}} \) and \(w _{i}^{\text {sm-p}} \), respectively. Note that these variables are the only ones having nonzero objective coefficients in case that the solution values of all slack variables are zero. The linear program \(\text {LP}_{\text {base}}(S)\) together with constraints (53)–(58) we denote as smoothing linear program \(\text {LP}_{\text {sm}}(S)\) in the following.
6 Iterative velocity adjustment procedure
A main drawback of the linear model for the transient gas flow through pipelines used in our MILP formulation is the fixation of the absolute velocity in the friction term of the Momentum Equation, see Sect. 4.8. If the mass flows or pressures at the endnodes of a pipeline change significantly over time, we might under- or overestimate the pressure loss, what in turn may lead to inaccurate decisions. Thus, we introduce an iterative velocity adjustment procedure (IVAP) which we later on combine with the hierarchical optimization model, see Sect. 7. Inspired by the successful application of sequential linear programming in the context of gas network control problems, see Rueda et al. (2019) for an example, its goal is to determine a solution, which is feasible for the tri-level MILP L featuring non-linear constraints (M) instead of the linear model (11) for the Momentum Equations.
The main idea is the following: Given a solution S for the hierarchical MILP model, we first of all consider the induced linear program \(\text {LP}_{\text {base}}(S)\) as defined in Sect. 5.2. Next, we retrieve the gas velocities at the endnodes of all pipelines w.r.t. to S and update the friction terms of the linear model for the Momentum Equations accordingly. The goal of this LP shall be to determine a solution such that the pressure and flow variables corresponding to constraints (11) stay as close as possible to the variables values in S. In case that there exists a solution such that all these variable values are equal, we have determined a solution satisfying constraints (M), as the gas velocity depends only on the mass flow and the pressure in our model.
Thus, we need to introduce additional variables and constraints to measure deviations and penalize them. Therefore, let us denote the set of vertices incident to a pipeline by \({\mathcal {V}}^{\text {pi}} {:=} \bigcup \nolimits _{(\ell ,r) \in {\mathcal {A}} ^{\text {pi}}} \{\ell ,r\}\). Furthermore, given a solution S for the tri-level MILP formulation, let \(p _{v,t} ^S\) \(q _{a,\ell ,t} ^S\), and \(q _{a,r,t} ^S\) denote the values of the corresponding variables in the following.
First of all, we add variables \(\delta _{v,t}^{p +},\delta _{v,t}^{p-} \in {\mathbb {R}}_{\ge 0} \) for each node \(v \in {\mathcal {V}}^{\text {pi}} \) and each time step \(t \in {\mathcal {T}} \). Furthermore, we add one additional variable \({\overline{\delta }}^{p}\) and constraints
Constraint (59) measures the difference of the pressure value of node v and time step t to the corresponding solution value in S using variables \(\delta _{v,t}^{p +}\) and \(\delta _{v,t}^{p-}\). The maximum difference for any node and any time step is determined by constraint (60) and is equal to variable \({\overline{\delta }}^{p}\), since it has non-zero objective coefficient \(w ^{\text {sm-p}} \in {\mathbb {R}}_{\ge 0} \).
Second, for each pipeline \(a = (\ell ,r) \in {\mathcal {A}} ^{\text {pi}} \) and each time step \(t \in {\mathcal {T}} \) we add four continuous variables \(\delta _{a,\ell ,t}^{q +}, \delta _{a,\ell ,t}^{q-}, \delta _{a,r,t}^{q +}, \delta _{a,r,t}^{q-} \in {\mathbb {R}}_{\ge 0} \) as well as one additional variable \({\overline{\delta }}^{q}\in {\mathbb {R}}_{\ge 0} \) and introduce constraints
Here, constraints (61) and (62) determine the difference between the in- or outflow at the endnodes \(\ell \) and r of pipeline a and the corresponding solution values in S. The maximum difference for any pipe, endnode and time step is determined by constraints (63) and (64) and is equal to \({\overline{\delta }}^{q}\)since it has non-zero objective coefficient \(w ^{\text {sm-q}} \in {\mathbb {R}}_{\ge 0} \).
The iterative velocity adjustment procedure (IVAP) can now be stated as shown in Algorithm 3: Given a feasible solution \(S_0\) for the initial hierarchical MILP model L, we first of all determine the gas velocities for all \(v \in {\mathcal {V}}^{\text {pi}} \) and all time steps w.r.t. \(S_0\). Next, we iteratively repeat the following procedure: In iteration i, we are going to use the absolute value of the average gas velocities from the last \(\min \{i,k\}\) solutions, i.e., \(|v_i^*|\), in constraints (11). Thus, we obtain linear program \(\text {LP}_{\text {iv}}^i\) as \(\text {LP}_{\text {base}}(S_{i-1})\) using \(|v_i^*|\) in constraints (11) together with constraints (59)–(64). If \(\text {LP}_{\text {iv}}^i\) is infeasible, the procedure is terminated and UNSUCCESSFUL is returned. Otherwise we retrieve the gas velocities \(v_i\) from an optimal solution \(S_i\). If \(|v_i|\) and \(|v_i^*|\) differ by less than \(\varepsilon \) for all pipelines, nodes and time steps, \(S_i\) is returned as result.
7 Algorithmic framework
The complete algorithmic framework for the optimal operation of transient gas transport networks, which combines AATM and the solution smoothing from Sect. 5 with the IVAP from Sect. 6 is shown in Algorithm 4.
The following procedure is iteratively repeated: In iteration i, we solve the tri-level MILP formulation L with the AATM. If the model is infeasible, UNSUCCESSFUL is returned. Otherwise, we apply the smoothing routine to its solution \(S_i\) and start the IVAP with the resulting \(S^{sm}_i\). If the IVAP terminates with a feasible solution, this is a feasible solution for L with constraints (M) instead of (11) and Algorithm 4 is terminated. Otherwise, we retrieve the velocities \(v^{\text {update}}_i\) from the last IVAP iteration and derive a new tri-level MILP formulation L where constraints (11) are updated using the corresponding absolute values in the friction terms. Further, we add a no-good-cut w.r.t. the simple state variable values in \(S_i\) to L, i.e.,
where \(x_{s,t} ^{S_i}\) denotes the corresponding variable value in \(S_i\), \(k+1\) is the number of time steps, and m the number of network stations. The algorithm terminates with UNSUCCESSFUL if no solution within \(\varDelta \) iterations is found.
The reason to start the IVAP with the smoothed solution \(S^{sm}\) instead of S is that the differences w.r.t. the gas velocities at the nodes for consecutive time steps are intuitively smaller. This may not only accelerate the procedure, but also lead to a higher success rate of the IVAP. Furthermore, the solution produced by the IVAP may maintain some degree of smoothness of the initial \(S^{sm}\).
If the IVAP fails, we use the velocities from its last iteration in the friction terms of the Momentum Equations of model L in the next iteration, since it may be closer to the velocities of some feasible solution for L with non-linear equations (M). By adding the no-good-cut ensures, that the same solution w.r.t. the simple states does not occur again. We do this, since we were not able to find a feasible solution for the non-linear pipe equations using this solution.
Finally, it is important to note, that the solutions produced here are not necessarily optimal for L with the non-linear equations (M). However, if no slack is needed and the objective value is zero, the solution is actually optimal.
8 Computational experiments
In this section we present computational experiments, which were conducted in order to test whether the framework presented in Algorithm 4 is suitable to make important transient global control decisions for the first stage of the KOMPASS algorithm or not. First, we describe our set of test instances, followed by the computational setup. We conclude with an analysis of the results.
8.1 Instances
For our computations we used 333 instances provided by our project partner OGE, which are based on a real-world network and the corresponding historically measured data. Thus, the initial state, the boundary values, and the source pressures represent feasible network states. Non-technical control decisions, which were undertaken by the dispatchers during the considered time horizon, are already included in this data. Hence, if the model would perfectly capture reality, we expect it to find feasible solutions without using any slack.
While the overall composition of the network is depicted in Table 1, the properties of the 7 network stations contained in the network are shown in Table 2. For our experiments, we considered two different data sets. Dataset 1 consists of 168 instances in 30 minute intervals starting at noon of a virtual day 1 and ending of 23:30pm of virtual day 4. The other set consists of 165 instances starting at midnight of a virtual day 6 and ending at 10:00am on virtual day 9. Note that the time difference of 30 minutes between two successive instances does not posses any meaning and is due to the data creation process at OGE. Furthermore, we considered a granularity of \(4 \cdot 15\) minutes and \(11 \cdot 60\) minutes, i.e., we have \(n=15\) time steps and cover a time horizon of 12 h for each single instance. The reason for the first four time steps having short length is to have more exact recommendations w.r.t. the time of their execution in the near future. Additionally, this goes hand in hand with one of our our long-term goals, namely to have an overall run time of KOMPASS of 15 minutes.
The cost parameters regarding the third level objective were set to \(w ^a = 5.0\) for all \(a \in {\mathcal {A}} ^{\text {ar}} \) and for each simple state \(s \in {\mathcal {S}} \) an individual cost parameter \(w ^s \in \{0,\ldots ,200\}\) for a change into it was fixed according to expert opinions. Recall, that the cost for both types of slack variables is equal to 1.0, see Sect. 4.13. For the solution smoothing we set \(w _{v}^{\text {sm-q}} {:=} 15\), \(w _{v}^{\text {sm-p}} {:=} 150\), \(w _{i}^{\text {sm-q}} {:=} |{\mathcal {V}} ^{\text {fn}} _{i} |\), and \(w _{i}^{\text {sm-p}} {:=} 10.0 |{\mathcal {V}} ^{\text {fn}} _{i} |\) for each \(v \in {\mathcal {V}} ^{\text {fn}} _{i} \) and network station \(G_{i} \), while for the IVAP we used \(w ^{\text {sm-p}} {:=} 10^4\) and \(w ^{\text {sm-q}} {:=} 10^3\), respectively.
8.2 Computational setup
We performed our computations on a cluster of machines composed of two Intel Xeon Gold 5122 running at 3.60 GHz, which provide in total 8 cores and 96 GB of RAM. As solver for the underlying MILP problems we used Gurobi in version 9.0.2 (Gurobi Optimization 2020), which was accessed via the native C interface.
Since the corresponding MILP formulations turned out to be numerically challenging, we set the NumericFocus parameter to the maximum value, as well as IntFeasTol and MIPGap to \(10^{-6}\), and used the standard settings of Gurobi otherwise. Additionally, we fixed the absolute velocities in the friction-based pressure difference term of the Momentum Equation (M), i.e., in the third summand, to the maximum of the absolute gas velocity in the initial time step at the corresponding node and \(v^{\mathrm{min}} {:=} 0.1 \frac{m}{s}\). The introduction of the threshold \(v^{\mathrm{min}}\) is necessary in order to control the magnitudes of the constraint’s coefficients and thereby to avoid numerical instabilities. However, this threshold is decreased to \(10^{-3}\frac{m}{s}\) for the IVAP, where we used \(k {:=} 3\) and \(\varepsilon {:=} 10^{-2}\frac{m}{s}\) as input parameters.
Furthermore, to each test instance we applied two versions of the AATM, i.e., Algorithm 1, which we used within the algorithmic framework, i.e., Algorithm 4. First, a pure approach, as described in the pseudocode (PURE). And second, a version where both the RHH and the MCF heuristic were run before every MILP of the overall solving process (HEUR). In both cases, we set a cumulative time limit of 3600 s for all MILPs that are solved within all calls of the AATM. For HEUR, we additionally imposed individual sub time limits of 300 s for each single MILP solve performed by RHH or MCF. For the smoothing LPs and the IVAP no time limit was imposed.
8.3 Results
Detailed computational results for all our test runs can be found in Appendix A. First, we compare the performance of the PURE and the HEUR approach w.r.t. the determined solutions and the cumulative run times of all AATM calls, see Fig. 4. Afterwards we discuss the run times of the whole framework, see Algorithm 4 and Fig. 5, and conclude with an analysis of the solutions and the objective values.
The y-axis represent time in seconds on a logarithmic scale. For each problem instance in data set 1 (left) and data set 2 (right), there exists a blue and an orange dot representing the cumulative running times of the PURE and the HEUR approach spent in the AATM. Note that dots can overlap if too many are present for the same value range
Considering the first data set, both approaches determined solutions with the same objective values in rather short amounts of time for all instances. Although it seems that PURE outperforms HEUR here, this is only true for the instances having zero objective value. Here, HEUR is slower as it has to run our heuristics, which take up the biggest shares of the running times. Additionally, feasible solutions having zero objective seem to be quickly found by Gurobi, too. However, on fourteen of the sixteen instances with non-zero objective HEUR is faster.
In the second data set, for all but one instance, namely 7-2230, a feasible solution satisfying the \(\varepsilon \)-criterion was found by the HEUR approach. In contrast, the PURE approach failed to determine any feasible solution during the first run of AATM for instances 6-0500, 6-0730, 6-0830, and 6-1400, and could not produce a solution satisfying the \(\varepsilon \)-criterion for instances 6-0800, 7-1630, and, analogous to HEUR, 7-2230. HEUR outperformed PURE w.r.t. the determined solution values, as it always terminated with not greater objective value, except for instance 6-1230. Although the PURE approach again solved instances with solution value zero faster, the advantage of HEUR on instances with non-zero objective becomes inevitably clear considering the upper part of the plot.
The y-axis represent time in seconds on a logarithmic scale. For each problem instance in data set 1 (left) and data set 2 (right), there exists a blue and an orange dot representing the running times of the PURE and the HEUR approach including smoothing and IVAP, respectively. Note that dots may overlap
If we consider the running times of the complete algorithmic approaches, which are shown in Fig. 5, the advantage of PURE becomes less significant as the IVAP takes up a big share of the complete run times for instances with zero objective. Additionally, we see that HEUR is the more robust approach, not only w.r.t. determining feasible solutions as discussed above, but also w.r.t. the run times, in particular when considering the second data set. For all but instances 6-1230 and 7-2230, HEUR terminated within the time limit of 3600 s even if smoothing and IVAP times are included.
Only one instance needed the inflow slack at boundary nodes: Instance 6-1530 had small amount of slack used in the second time step. A deeper analysis showed, that this slack was needed close to a network station, suggesting that a deeper investigation regarding the corresponding station and the input data is necessary. Furthermore, for a majority of the instances a solution was found within the first iteration of the framework. For HEUR, for example, this is the case for 315 instances. Additionally, the IVAP performed around 50–70 iterations for most of the instances.
Considering the solution value curves shown in Fig. 6 and analyzing the corresponding measures, we can see that the recommended technical control measures are stable over consecutive runs, which is a desired behaviour for our decision support system KOMPASS. Most importantly, the experts consider them reasonable. However, for the outliers, in particular instances 4-0700, 6-0530, 6-0600, and 6-2000 a deeper investigation of the model and the corresponding input data has to be conducted.
While the y-axis shows the third level’s objective value of the solution, the x-axis represents the problem instances of data set 1 (left) and data set 2 (right) in chronological order. The curves are a linear interpolation of the values of the best solution w.r.t. the third level that was found by any of the two approaches, i.e., all HEUR solutions except for 6-1230. No value for instance 7-2230 is shown as no solution satisfying the \(\varepsilon \)-criterion was found by any of the two approaches
9 Conclusion and outlook
In this paper we presented an algorithmic framework consisting of a tri-level MILP model together with a sequential linear programming inspired post-processing procedure for the optimization of the control of transient gas transport networks. Complex pipeline intersection areas are replaced by so-called network stations, i.e., simplified hand-tailored graph representations modelling the technical control capabilities. The proposed framework represents the first stage of the two-stage approach implemented in KOMPASS, a decision support system giving technical and non-technical control recommendations to the dispatchers. Its goal is to make important transient global control decisions, i.e., to determine the directions of the flow and where to compress the gas. The results are then verified in the second stage using the highly detailed model described in Hennings et al. (2019). Using a linear model for the transient gas flow in pipelines and further approximations, especially regarding the compressor stations, we derive a tri-level MILP formulation, which can be solved using a sequence of single-level MILP models. For these single-level MILPs, we developed two heuristics, which determine initial solutions of good quality in rather short amounts of time and have a positive impact on the overall run time. Additionally, we introduced an iterative velocity adjustment procedure IVAP inspired by sequential linear programming in order to eliminate the drawbacks of the linear model for the pipe equations. Our computational experiments using historic flow and pressure values suggest that the presented framework represents a valuable basis for further development and possible extensions of KOMPASS.
Regarding the model itself, the analysis of our computational results showed that a deeper investigation of the solutions w.r.t. the hand-tailored network stations has to follow. In particular, an automated process to create the network station models given the topologies of the intersection areas and their corresponding operation modes is necessary, as it would not only improve the robustness and accuracy of our framework, but also lead to a better maintainability of KOMPASS. This is because currently every time the topology at a junction changes, the network stations have to be revised and adapted by hand.
Additionally, there exists further potential regarding the modelling of the compressing arcs. Besides an improvement of the approximation of the power bound, it would be beneficial to dynamically adapt the compression ratio in constraint (40). Currently, we are working on including this feature in the IVAP.
Another major point for future research is a theoretical analysis of the IVAP. Although it works well in practice, it is our goal to derive some certificate of convergence. As a first step, we are currently working on a proof for passive networks, i.e., networks consisting of pipelines only, and for solutions of the tri-level MILP that fulfill certain criteria. Furthermore, we investigate if and how smoothing may be integrated.
From the industrial application point of view, we currently aim at applying the model to larger parts of the network, which include more intersection areas. Of course, it may also be beneficial to extend the formulation by modelling parallel and sequential setting choices for compressor stations or to improve the overall solution quality by adding additional features like ramp up times of compressor units. However, our current priority is to decrease the run time in order to be able to scale the presented approach. Hence, we are for example investigating possibilities to speed up the solving process of the single-level MILPs in the AATM by developing new heuristics and the introduction of valid cutting planes. Additionally, we are looking into speeding up the IVAP by developing a warm start heuristic for successive LP runs.
The whole algorithmic framework and the IVAP in particular may also benefit from a new feature that we are currently integrating: Using the solutions of chronologically earlier instances. Assume you have solved some problem instance, consider the next instance based on the network state 30 minutes later. First of all, we are working on a new primal heuristic based on the solution values of the binary variables corresponding to the jointly covered time horizon. And second, we are experimenting with using the absolute velocities of the solution directly in the Momentum Eq. (11) of the initial tri-level MILP model for the new instance. The rationale here is, that those values probably serve as a better and more realistic starting point.
Finally, for our computational experiments we used instances that are based on historic supplies, demands and pressure values. When forecasted values are used, the presented model may need to apply non-technical measures more often. This can lead to a significant increase in the run time. Thus, this is an issue that we have to address in the near future, too.
References
Achterberg T (2007) Constraint Integer Programming. Doctoral thesis, Technische Universität Berlin, Fakultät II - Mathematik und Naturwissenschaften, Berlin. https://doi.org/10.14279/depositonce-1634. URL http://dx.doi.org/10.14279/depositonce-1634
Addis B, Carello G, Grosso A, Tànfani E (2016) Operating room scheduling and rescheduling: a rolling horizon approach. Flex Serv Manuf J 28(1–2):206–232
Alderson DL, Brown GG, Carlyle WM, Wood RK (2011) Solving Defender-Attacker-Defender Models for Infrastructure Defense. Operations Research, Computing, and Homeland Defense, INFORMS pp 28–49
Brotcorne L, Labbé M, Marcotte P, Savard G (2001) A bilevel model for toll optimization on a multicommodity transportation network. Transp Sci 35(4):345–358
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/s11081-018-9414-5
Colson B, Marcotte P, Savard G (2007) An overview of bilevel optimization. Ann Oper Res 153(1):235–256
Dempe S (2002) Foundations of bilevel programming. Springer Science & Business Media, Berlin
Dempe S, Kalashnikov V, Ríos-Mercado RZ (2005) Discrete bilevel programming: Application to a natural gas cash-out problem. Eur J Oper Res 166(2):469–488
Domschke P, Geißler B, Kolb O, Lang J, Martin A, Morsi A (2011) Combination of Nonlinear and Linear Optimization of Transient Gas Networks. Inform J Comput 23(4):605–617
Federal Ministry for Economic Affairs and Energy: Still indispensable for a reliable energy supply. https://www.bmwi.de/Redaktion/EN/Dossier/conventional-energy-sources.html (2019). Accessed: 2019-10-31
Forschungscampus MODAL. http://forschungscampus-modal.de/. Accessed: 2019-10-31
Fügenschuh A, Geißler B, Gollmer R, Morsi A, Pfetsch ME, Rövekamp J, Schmidt M, Spreckelsen K, Steinbach MC. Physical and technical fundamentals of gas networks. In: Koch et al. [26], (2015)
Gao Z, Wu J, Sun H (2005) Solution algorithm for the bi-level discrete network design problem. Transp Res Part B: Methodol 39(6):479–495
Garcia-Herreros P, Zhang L, Misra P, Arslan E, Mehta S, Grossmann IE (2016) Mixed-integer bilevel optimization for capacity planning with rational markets. Comput Chem Eng 86:33–47
Grimm V, Schewe L, Schmidt M, Zöttl G (2019) A multilevel model of the European entry-exit gas market. Math Methods Oper Res 89(2):223–255
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/s10589-017-9970-1
Gurobi Optimization L (2020) Gurobi optimizer reference manual, version 9.0.2. http://www.gurobi.com
Hennings F (2018) Benefits and Limitations of Simplified Transient Gas Flow Formulations. In: Operations research proceedings 2017, pp 231–237. Springer
Hennings F, Anderson L, Hoppmann K, Turner M, Koch T (2020) Controlling transient gas flow in real-world pipeline intersection areas. Optim Eng, Springer, pp 1–48
Hoppmann K, Schwarz R (2018) Finding maximum minimum cost flows to evaluate gas network capacities. Oper Res Proc 2017:339–346
Jiang R, Wang J, Guan Y (2011) Robust unit commitment with wind power and pumped storage hydro. IEEE Trans Power Syst 27(2):800–810
Kalashnikov VV, Dempe S, Pérez-Valdés GA, Kalashnykova NI, Camacho-Vallejo JF (2015) Bilevel Programming and Applications. Mathematical Problems in Engineering 2015
Kalashnikov VV, Ríos-Mercado RZ (2002) Solving a natural gas cash-out bilevel program by a penalty function method. Technical Report PISIS-2002-03, Graduate Program in Systems Engineering, UANL, San Nicolás de los Garza, México
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. Numer Algorithms 53(2–3):293–307
Lu J, Han J, Hu Y, Zhang G (2016) Multilevel decision-making: a survey. Inf Sci 346:463–487
Mak TW, Van Hentenryck P, Zlotnik A, Hijazi H, Bent R (2016) Efficient dynamic compressor optimization in natural gas transmission systems. In: American control conference (ACC), 2016, pp 7484–7491. IEEE. https://doi.org/10.1109/ACC.2016.7526855
Migdalas A, Pardalos PM, Värbrand P (2013) Multilevel optimization: algorithms and applications, vol 20. Springer Science & Business Media, Berlin
Moritz S (2007) A mixed integer approach for the transient case of gas network optimization. Ph.D. thesis, Technische Universität Darmstadt, Darmstadt
Network data of Open Grid Europe GmbH. https://www.open-grid-europe.com/cps/rde/oge-internet/hs.xsl/Strukturdaten-gemass-27-Abs-2-GasNEV-654.htm?rdeLocaleAttr=de. Accessed: 2019-10-31
Nielsen LK, Kroon L, Maróti G (2012) A rolling horizon approach for disruption management of railway rolling stock. Eur J Oper Res 220(2):496–509
Nikuradse J (1950) Laws of Flow in Rough Pipes. The National Advisory Committee for Aeronautics (NACA), Washington
Open Grid Europe GmbH. https://www.open-grid-europe.com/. Accessed: 2019-10-31
Osiadacz AJ (1996) Different transient flow models - limitations, advantages, and disadvantages. Paper presented at the PSIG Annual Meeting, San Francisco, California
Petkovic M, Chen Y, Gamrath I, Gotzes U, Hadjidimitriou NS, Zittel J, Xu X, Koch T (2019) A Hybrid Approach for High Precision Prediction of Gas Flows. Technical Report pp 19–26, ZIB, Takustr. 7, 14195 Berlin
Rövekamp J (2015) Background on gas market regulation. In: Koch et al. [26]
Rueda ÁMG, Díaz JG, de Córdoba MPF (2019) A twist on slp algorithms for nlp and minlp problems: an application to gas transmission networks. Optim Eng 20(2):349–395
Saleh JM (2002) Fluid Flow Handbook. McGraw-Hill Professional, USA
Samà M, D’Ariano A, Pacciarelli D (2013) Rolling horizon approach for aircraft scheduling in the terminal control area of busy airports. Proc-Soc Behav Sci 80:531–552
Steringa JJ, Hoogwerf M, Dijkhuis H (2015) A Systematic Approach to Transmission Stress Tests in Entry-Exit Systems. Paper presented at the PSIG Annual Meeting, New Orleans, Louisiana
Winegarden W (2019) Low- and no-carbon future starts with natural gas. https://www.washingtonpost.com/brand-studio/wp/2019/02/15/low-and-no-carbon-future-starts-with-natural-gas/. Accessed: 2019-10-31
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), pp 4563–4570. IEEE. https://doi.org/10.1109/CDC.2015.7402932
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 numbers 05M14ZAM, 05M20ZBM).
Funding
Open Access funding enabled and organized by Projekt DEAL.
Author information
Authors and Affiliations
Corresponding author
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Appendix
Appendix
The columns of the following tables contain the results of the computational experiments described and discussed in Sect. 8. In the following, the solution values are stated w.r.t. the objective function of the third level (52), i.e., the cost of slack is not shown. All running times are given in seconds.
The 1st column contains the instance name. It consists of the virtual day together with the time in hours and minutes of the corresponding initial state, i.e., 2-0400 is the instance having the initial state from 4am of virtual day 2. For cursively written instances both approaches featured some non-zero solution values for the inflow slack variables. Otherwise all slack variables were zero for the final solutions.
The 2nd and 3rd column show the value of the best solution found and the gap for the final MILP solve within AATM, respectively. In case these columns are empty, the approach did not find a feasible solution in its last iteration.
The 3rd column denotes the total run time spent in the algorithmic approach AATM, i.e., Algorithm 1 for the PURE variant, while the 4th column denotes the run time of the complete algorithmic framework, i.e., it additionally includes solution smoothing and the runs of the iterative velocity adjustment procedure (IVAP). If the time limit of 3600 s for the AATM runs was hit and the 4th column also lists 3600 s as result, this means that no feasible solution was found by AATM within the time limit in the first iteration of Algorithm 4 since no smoothing or velocity adjustment took place. Otherwise, either the best solution found was smoothed and processed by the IVAP or at least one no-good-cut must have been added, i.e., Algorithm 4 performed at least two iterations.
Additionally, the 6th column shows the number of total IVAP iterations, while the 7th column shows the number of added no-good-cuts which is equal to the number of iterations of Algorithm 4. The 8th up to the 13th column contain the corresponding entities for the HEUR approach. Note that the run times of the heuristics are included in the 10th column.
INST | Val | Gap | AATM | PURE | SM-IT | NG | Val | Gap | AATM | HEUR | SM-IT | NG |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
1-1200 | 0 | 0 | 10.99 | 225.98 | 56 | 0 | 0 | 0 | 23.1 | 326.21 | 54 | 0 |
1-1230 | 0 | 0 | 13.36 | 236.96 | 44 | 0 | 0 | 0 | 23.87 | 238.17 | 43 | 0 |
1-1300 | 0 | 0 | 12.62 | 211.04 | 48 | 0 | 0 | 0 | 52.25 | 229.61 | 45 | 1 |
1-1330 | 0 | 0 | 14.89 | 263.45 | 44 | 0 | 0 | 0 | 19.09 | 185.72 | 49 | 0 |
1-1400 | 0 | 0 | 12.43 | 377.97 | 47 | 0 | 0 | 0 | 21.15 | 278.8 | 50 | 0 |
1-1430 | 0 | 0 | 12.77 | 217.51 | 59 | 0 | 0 | 0 | 24.92 | 221.21 | 52 | 0 |
1-1500 | 0 | 0 | 12.65 | 339.93 | 46 | 0 | 0 | 0 | 27.23 | 490.1 | 60 | 0 |
1-1530 | 0 | 0 | 12.28 | 377.6 | 42 | 0 | 0 | 0 | 18.74 | 216.09 | 43 | 0 |
1-1600 | 0 | 0 | 11.12 | 212.29 | 45 | 0 | 0 | 0 | 19.79 | 177.99 | 43 | 0 |
1-1630 | 0 | 0 | 12.4 | 197.8 | 48 | 0 | 0 | 0 | 19.26 | 175.84 | 46 | 0 |
1-1700 | 0 | 0 | 10.99 | 181.47 | 43 | 0 | 0 | 0 | 19.51 | 179.2 | 46 | 0 |
1-1730 | 0 | 0 | 11.67 | 248.37 | 53 | 0 | 0 | 0 | 24.5 | 195.33 | 45 | 0 |
1-1800 | 0 | 0 | 67.32 | 275.53 | 54 | 0 | 0 | 0 | 23.82 | 211.6 | 46 | 0 |
1-1830 | 0 | 0 | 11.44 | 179.18 | 45 | 0 | 0 | 0 | 25.5 | 208.27 | 47 | 0 |
1-1900 | 0 | 0 | 12.13 | 176.09 | 46 | 0 | 0 | 0 | 24.89 | 317.58 | 50 | 0 |
1-1930 | 0 | 0 | 19.21 | 202.26 | 49 | 0 | 0 | 0 | 39.8 | 223.73 | 49 | 0 |
1-2000 | 0 | 0 | 12.87 | 197.46 | 47 | 0 | 0 | 0 | 22.25 | 226.12 | 49 | 0 |
1-2030 | 0 | 0 | 12.98 | 280.37 | 45 | 0 | 0 | 0 | 22.76 | 193.94 | 49 | 0 |
1-2100 | 0 | 0 | 11.35 | 254.89 | 61 | 0 | 0 | 0 | 24.33 | 202.53 | 46 | 0 |
1-2130 | 0 | 0 | 13 | 187.86 | 50 | 0 | 0 | 0 | 21.86 | 190.06 | 46 | 0 |
1-2200 | 0 | 0 | 11.65 | 200.2 | 51 | 0 | 0 | 0 | 19.12 | 165.17 | 43 | 0 |
1-2230 | 0 | 0 | 11.55 | 222.28 | 44 | 0 | 0 | 0 | 19.66 | 282.64 | 50 | 0 |
1-2300 | 0 | 0 | 13.95 | 200.91 | 45 | 0 | 0 | 0 | 22.3 | 218.06 | 49 | 0 |
INST | Val | Gap | AATM | PURE | SM-IT | NG | Val | Gap | AATM | HEUR | SM-IT | NG |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
1-2330 | 0 | 0 | 11.97 | 200.11 | 50 | 0 | 0 | 0 | 22.39 | 216.22 | 47 | 0 |
2-0000 | 0 | 0 | 11.78 | 210.78 | 50 | 0 | 0 | 0 | 24.77 | 326.06 | 55 | 0 |
2-0030 | 0 | 0 | 13.39 | 204.14 | 51 | 0 | 0 | 0 | 26.51 | 218.03 | 46 | 0 |
2-0100 | 0 | 0 | 12.53 | 188.86 | 44 | 0 | 0 | 0 | 24.43 | 365.18 | 49 | 0 |
2-0130 | 0 | 0 | 12.7 | 235.53 | 53 | 0 | 0 | 0 | 32.54 | 253.63 | 53 | 0 |
2-0200 | 0 | 0 | 12.58 | 286.92 | 49 | 0 | 0 | 0 | 19.5 | 163.44 | 45 | 0 |
2-0230 | 0 | 0 | 12.93 | 241.62 | 49 | 0 | 0 | 0 | 33.12 | 260.76 | 49 | 0 |
2-0300 | 0 | 0 | 12.3 | 310.6 | 53 | 0 | 0 | 0 | 20.59 | 381.86 | 54 | 0 |
2-0330 | 0 | 0 | 12.22 | 267.06 | 50 | 0 | 0 | 0 | 19.69 | 334.53 | 62 | 0 |
2-0400 | 0 | 0 | 11.37 | 200.32 | 46 | 0 | 0 | 0 | 20.11 | 202.44 | 51 | 0 |
2-0430 | 0 | 0 | 12.96 | 248.93 | 51 | 0 | 0 | 0 | 20.58 | 283.96 | 47 | 0 |
2-0500 | 0 | 0 | 14.93 | 232.7 | 57 | 0 | 0 | 0 | 26.46 | 227.23 | 54 | 0 |
2-0530 | 0 | 0 | 12.26 | 218.2 | 48 | 0 | 0 | 0 | 26.85 | 193.18 | 43 | 0 |
2-0600 | 0 | 0 | 12.31 | 208.37 | 52 | 0 | 0 | 0 | 23.48 | 201.94 | 47 | 0 |
2-0630 | 0 | 0 | 13.83 | 233.46 | 46 | 0 | 0 | 0 | 26.32 | 211.75 | 45 | 0 |
2-0700 | 0 | 0 | 12.1 | 219.65 | 46 | 0 | 0 | 0 | 26.12 | 194.7 | 44 | 0 |
2-0730 | 0 | 0 | 12.88 | 278.35 | 48 | 0 | 0 | 0 | 26.61 | 221.86 | 46 | 0 |
2-0800 | 0 | 0 | 12.02 | 228.06 | 41 | 0 | 0 | 0 | 28.69 | 192.77 | 43 | 0 |
2-0830 | 0 | 0 | 11.83 | 232.01 | 41 | 0 | 0 | 0 | 26.52 | 193.45 | 45 | 0 |
2-0900 | 0 | 0 | 15.12 | 201.12 | 49 | 0 | 0 | 0 | 31.95 | 209.63 | 51 | 0 |
2-0930 | 0 | 0 | 14.07 | 257.42 | 44 | 0 | 0 | 0 | 33.35 | 179.89 | 42 | 0 |
2-1000 | 0 | 0 | 13.42 | 215.57 | 46 | 0 | 0 | 0 | 28.73 | 189.12 | 46 | 0 |
2-1030 | 0 | 0 | 11.47 | 205.77 | 43 | 0 | 0 | 0 | 27.12 | 180.32 | 43 | 0 |
2-1100 | 0 | 0 | 11.92 | 189.32 | 44 | 0 | 0 | 0 | 23.46 | 186.7 | 45 | 0 |
2-1130 | 0 | 0 | 12.78 | 230.85 | 53 | 0 | 0 | 0 | 316.24 | 483.59 | 44 | 0 |
2-1200 | 0 | 0 | 12.59 | 287.21 | 45 | 0 | 0 | 0 | 632.67 | 847.09 | 52 | 1 |
2-1230 | 0 | 0 | 12.35 | 177.1 | 46 | 0 | 0 | 0 | 313.37 | 493 | 46 | 0 |
2-1300 | 0 | 0 | 11.89 | 198.78 | 47 | 0 | 0 | 0 | 313.32 | 526.14 | 48 | 0 |
2-1330 | 0 | 0 | 11.75 | 264.54 | 54 | 0 | 0 | 0 | 283.86 | 436.63 | 43 | 0 |
2-1400 | 0 | 0 | 11.91 | 233.64 | 53 | 0 | 0 | 0 | 313.09 | 473.44 | 48 | 0 |
2-1430 | 0 | 0 | 11.07 | 185.68 | 45 | 0 | 0 | 0 | 313.1 | 470.03 | 42 | 0 |
2-1500 | 0 | 0 | 11.18 | 203.62 | 52 | 0 | 0 | 0 | 310.11 | 470.53 | 53 | 0 |
2-1530 | 0 | 0 | 10.49 | 175.66 | 46 | 0 | 0 | 0 | 64.13 | 217.83 | 40 | 0 |
2-1600 | 0 | 0 | 24.97 | 212.57 | 55 | 0 | 0 | 0 | 312.92 | 525.89 | 58 | 0 |
2-1630 | 0 | 0 | 10.9 | 218.19 | 58 | 0 | 0 | 0 | 182.26 | 345.83 | 46 | 0 |
2-1700 | 0 | 0 | 13.26 | 219.93 | 54 | 0 | 0 | 0 | 319.2 | 473.1 | 47 | 0 |
2-1730 | 0 | 0 | 11.93 | 208.27 | 48 | 0 | 0 | 0 | 43.21 | 203.75 | 45 | 0 |
2-1800 | 0 | 0 | 10.52 | 196.08 | 49 | 0 | 0 | 0 | 257.03 | 394.22 | 45 | 0 |
2-1830 | 0 | 0 | 10.94 | 182.13 | 46 | 0 | 0 | 0 | 92.07 | 241.15 | 41 | 0 |
2-1900 | 0 | 0 | 10.83 | 190.26 | 51 | 0 | 0 | 0 | 63.38 | 219.95 | 43 | 0 |
2-1930 | 0 | 0 | 10.66 | 181.32 | 52 | 0 | 0 | 0 | 55.8 | 400.2 | 57 | 0 |
2-2000 | 0 | 0 | 13.65 | 277.3 | 49 | 0 | 0 | 0 | 67.82 | 491.07 | 66 | 0 |
2-2030 | 0 | 0 | 10.68 | 371.18 | 50 | 0 | 0 | 0 | 76.76 | 397.1 | 48 | 0 |
2-2100 | 0 | 0 | 11.2 | 295.92 | 50 | 0 | 0 | 0 | 60.69 | 368.89 | 45 | 0 |
2-2130 | 0 | 0 | 11.96 | 284.45 | 49 | 0 | 0 | 0 | 50.44 | 451.92 | 54 | 0 |
2-2200 | 0 | 0 | 11.94 | 484.38 | 52 | 0 | 0 | 0 | 63.31 | 557.59 | 51 | 0 |
2-2230 | 0 | 0 | 12.21 | 384.72 | 61 | 0 | 0 | 0 | 82.51 | 458.42 | 62 | 0 |
2-2300 | 0 | 0 | 12.96 | 330.6 | 49 | 0 | 0 | 0 | 26.09 | 293.37 | 50 | 0 |
2-2330 | 0 | 0 | 11.78 | 327.02 | 52 | 0 | 0 | 0 | 23.94 | 444.13 | 54 | 0 |
3-0000 | 0 | 0 | 12.02 | 318.39 | 53 | 0 | 0 | 0 | 20.22 | 451.52 | 52 | 0 |
3-0030 | 0 | 0 | 11.97 | 313.32 | 50 | 0 | 0 | 0 | 22.69 | 439.95 | 54 | 0 |
3-0100 | 0 | 0 | 13.09 | 228.76 | 51 | 0 | 0 | 0 | 19.99 | 464.64 | 52 | 0 |
3-0130 | 0 | 0 | 13.88 | 437.81 | 55 | 0 | 0 | 0 | 19.51 | 286.08 | 82 | 0 |
3-0200 | 0 | 0 | 13.05 | 232.45 | 53 | 0 | 0 | 0 | 20.5 | 212.65 | 54 | 0 |
3-0230 | 0 | 0 | 13.96 | 232.06 | 56 | 0 | 0 | 0 | 19.01 | 191.04 | 65 | 0 |
3-0300 | 0 | 0 | 14.14 | 211.14 | 54 | 0 | 0 | 0 | 19.61 | 255.8 | 54 | 0 |
3-0330 | 0 | 0 | 13.94 | 195.78 | 50 | 0 | 0 | 0 | 19.6 | 219.28 | 54 | 0 |
3-0400 | 0 | 0 | 13.09 | 185.83 | 55 | 0 | 0 | 0 | 19.7 | 230.3 | 54 | 0 |
3-0430 | 0 | 0 | 13.7 | 169.11 | 53 | 0 | 0 | 0 | 19.77 | 222.12 | 52 | 0 |
3-0500 | 0 | 0 | 15.01 | 190.72 | 49 | 0 | 0 | 0 | 19.58 | 259.89 | 50 | 0 |
3-0530 | 0 | 0 | 12.95 | 206.87 | 55 | 0 | 0 | 0 | 21.28 | 231.68 | 51 | 0 |
INST | Val | Gap | AATM | PURE | SM-IT | NG | Val | Gap | AATM | HEUR | SM-IT | NG |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
3-0600 | 0 | 0 | 13.2 | 158.39 | 45 | 0 | 0 | 0 | 20.25 | 158.04 | 47 | 0 |
3-0630 | 0 | 0 | 12.12 | 184.94 | 49 | 0 | 0 | 0 | 19.8 | 205.16 | 53 | 0 |
3-0700 | 0 | 0 | 14.55 | 172.81 | 49 | 0 | 0 | 0 | 21.42 | 172.62 | 46 | 0 |
3-0730 | 0 | 0 | 15.86 | 198.43 | 41 | 0 | 0 | 0 | 20.24 | 404.88 | 48 | 0 |
3-0800 | 0 | 0 | 12.82 | 163.51 | 43 | 0 | 0 | 0 | 26.51 | 361.16 | 63 | 0 |
3-0830 | 0 | 0 | 12.84 | 164.7 | 43 | 0 | 0 | 0 | 26.49 | 361.26 | 63 | 0 |
3-0900 | 0 | 0 | 17.64 | 186.09 | 53 | 0 | 0 | 0 | 20.23 | 326.13 | 63 | 0 |
3-0930 | 0 | 0 | 14.28 | 151.34 | 44 | 0 | 0 | 0 | 20.51 | 179.17 | 43 | 0 |
3-1000 | 0 | 0 | 13.18 | 221.99 | 58 | 0 | 0 | 0 | 19.57 | 194.74 | 46 | 0 |
3-1030 | 0 | 0 | 14.62 | 161.69 | 44 | 0 | 0 | 0 | 20.46 | 173.49 | 46 | 0 |
3-1100 | 0 | 0 | 12.55 | 181.64 | 51 | 0 | 0 | 0 | 21.41 | 243.33 | 55 | 0 |
3-1130 | 0 | 0 | 14.33 | 167.84 | 44 | 0 | 0 | 0 | 18.96 | 161.62 | 45 | 0 |
3-1200 | 0 | 0 | 13.7 | 276.92 | 72 | 0 | 0 | 0 | 23.57 | 172.71 | 45 | 0 |
3-1230 | 0 | 0 | 15.3 | 168.5 | 43 | 0 | 0 | 0 | 26.16 | 196.97 | 45 | 0 |
3-1300 | 0 | 0 | 14.96 | 176.69 | 47 | 0 | 0 | 0 | 310.07 | 470.96 | 47 | 0 |
3-1330 | 0 | 0 | 14.1 | 165.41 | 42 | 0 | 0 | 0 | 75.81 | 240.29 | 43 | 0 |
3-1400 | 0 | 0 | 14.34 | 204.37 | 48 | 0 | 0 | 0 | 313.38 | 486.96 | 56 | 0 |
3-1430 | 0 | 0 | 15.73 | 175.93 | 44 | 0 | 0 | 0 | 117.06 | 281.22 | 44 | 0 |
3-1500 | 0 | 0 | 18.15 | 316.25 | 83 | 0 | 0 | 0 | 313.56 | 497.99 | 65 | 0 |
3-1530 | 0 | 0 | 13.78 | 208.88 | 46 | 0 | 0 | 0 | 270.71 | 458.73 | 45 | 0 |
3-1600 | 0 | 0 | 15.88 | 226.05 | 49 | 0 | 0 | 0 | 300.95 | 500.15 | 49 | 0 |
3-1630 | 0 | 0 | 17.87 | 243.05 | 43 | 0 | 0 | 0 | 267.39 | 521.95 | 43 | 0 |
3-1700 | 0 | 0 | 73.51 | 336.35 | 46 | 0 | 0 | 0 | 312.81 | 564.04 | 58 | 0 |
3-1730 | 0 | 0 | 13.12 | 251.13 | 45 | 0 | 0 | 0 | 24 | 260.87 | 44 | 0 |
3-1800 | 0 | 0 | 14.48 | 202.74 | 46 | 0 | 0 | 0 | 19.4 | 239.49 | 73 | 0 |
3-1830 | 0 | 0 | 14.36 | 267.46 | 49 | 0 | 0 | 0 | 23.89 | 231.77 | 43 | 0 |
3-1900 | 0 | 0 | 13.19 | 197.49 | 45 | 0 | 0 | 0 | 24.13 | 219.42 | 43 | 0 |
3-1930 | 0 | 0 | 14.53 | 249.2 | 44 | 0 | 0 | 0 | 43.42 | 250.42 | 41 | 0 |
3-2000 | 0 | 0 | 13.55 | 180.63 | 53 | 0 | 0 | 0 | 22.52 | 220.46 | 61 | 0 |
3-2030 | 0 | 0 | 15.28 | 213.26 | 46 | 0 | 0 | 0 | 22.71 | 207.46 | 45 | 0 |
3-2100 | 105 | 0 | 26.73 | 265.24 | 49 | 0 | 105 | 0 | 24.19 | 204.8 | 44 | 0 |
3-2130 | 105 | 0 | 22.45 | 397.37 | 56 | 0 | 105 | 0 | 24.42 | 239.2 | 49 | 0 |
3-2200 | 105 | 0 | 110.05 | 270.53 | 49 | 0 | 105 | 0 | 25.58 | 198.52 | 45 | 0 |
3-2230 | 105 | 0 | 30.26 | 196.31 | 43 | 0 | 105 | 0 | 27.95 | 216.4 | 48 | 0 |
3-2300 | 105 | 0 | 84.35 | 354.46 | 53 | 0 | 105 | 0 | 24.35 | 198.91 | 44 | 0 |
3-2330 | 105 | 0 | 141.72 | 416.5 | 53 | 0 | 105 | 0 | 27.05 | 302.85 | 52 | 0 |
4-0000 | 0 | 0 | 13.99 | 244.42 | 54 | 0 | 0 | 0 | 21.44 | 256.56 | 47 | 0 |
4-0030 | 0 | 0 | 13.06 | 244.35 | 50 | 0 | 0 | 0 | 20.68 | 292.13 | 60 | 0 |
4-0100 | 0 | 0 | 12.28 | 189.63 | 50 | 0 | 0 | 0 | 20.04 | 219.04 | 50 | 0 |
4-0130 | 0 | 0 | 14.28 | 202.48 | 50 | 0 | 0 | 0 | 25.64 | 230.92 | 48 | 0 |
4-0200 | 0 | 0 | 11.18 | 195.47 | 50 | 0 | 0 | 0 | 22.36 | 220.06 | 50 | 0 |
4-0230 | 0 | 0 | 14.1 | 227.4 | 49 | 0 | 0 | 0 | 20.66 | 199.52 | 48 | 0 |
4-0300 | 0 | 0 | 14.22 | 212.59 | 48 | 0 | 0 | 0 | 21.08 | 199.41 | 52 | 0 |
4-0330 | 0 | 0 | 12.26 | 187 | 48 | 0 | 0 | 0 | 20.09 | 230.91 | 49 | 0 |
4-0400 | 0 | 0 | 14.26 | 180.48 | 50 | 0 | 0 | 0 | 34.29 | 200.74 | 50 | 0 |
4-0430 | 0 | 0 | 13.96 | 182.07 | 48 | 0 | 0 | 0 | 19.69 | 225.56 | 48 | 0 |
4-0500 | 0 | 0 | 12.89 | 203.21 | 48 | 0 | 0 | 0 | 22.22 | 207.54 | 48 | 0 |
4-0530 | 0 | 0 | 13.81 | 190.98 | 61 | 0 | 0 | 0 | 19.58 | 233.13 | 64 | 0 |
4-0600 | 0 | 0 | 13.04 | 170.61 | 49 | 0 | 0 | 0 | 21.87 | 205.06 | 52 | 0 |
4-0630 | 0 | 0 | 28.59 | 218.78 | 46 | 1 | 0 | 0 | 21.38 | 180.5 | 49 | 0 |
4-0700 | 85 | 0 | 131.06 | 320.75 | 49 | 1 | 85 | 0 | 44 | 220.33 | 62 | 1 |
4-0730 | 0 | 0 | 17.4 | 184.51 | 57 | 0 | 0 | 0 | 20.14 | 173.61 | 50 | 0 |
4-0800 | 0 | 0 | 13.7 | 171.31 | 50 | 0 | 0 | 0 | 22.98 | 162.78 | 50 | 0 |
4-0830 | 0 | 0 | 12.97 | 177.93 | 51 | 0 | 0 | 0 | 22.59 | 191.5 | 51 | 0 |
4-0900 | 0 | 0 | 36.26 | 210.06 | 50 | 1 | 0 | 0 | 45.29 | 217.49 | 44 | 1 |
4-0930 | 0 | 0 | 14.34 | 176.17 | 46 | 0 | 0 | 0 | 20.62 | 189.69 | 47 | 0 |
4-1000 | 0 | 0 | 22.42 | 185.97 | 49 | 0 | 0 | 0 | 61.32 | 223.6 | 44 | 0 |
4-1030 | 0 | 0 | 13.8 | 157.16 | 48 | 0 | 0 | 0 | 20.26 | 166.68 | 42 | 0 |
4-1100 | 0 | 0 | 15.57 | 174.55 | 43 | 0 | 0 | 0 | 20.59 | 161.03 | 43 | 0 |
4-1130 | 0 | 0 | 13.33 | 174.29 | 48 | 0 | 0 | 0 | 21.14 | 164.76 | 42 | 0 |
INST | Val | Gap | AATM | PURE | SM-IT | NG | Val | Gap | AATM | HEUR | SM-IT | NG |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
4-1200 | 0 | 0 | 14.22 | 177.04 | 45 | 0 | 0 | 0 | 21.7 | 172.8 | 47 | 0 |
4-1230 | 0 | 0 | 14.48 | 237 | 66 | 0 | 0 | 0 | 22.18 | 183.07 | 46 | 0 |
4-1300 | 0 | 0 | 16.38 | 283.37 | 55 | 0 | 0 | 0 | 20.77 | 171.32 | 46 | 0 |
4-1330 | 0 | 0 | 14.28 | 305.25 | 58 | 0 | 0 | 0 | 21.41 | 241.34 | 44 | 0 |
4-1400 | 0 | 0 | 12.98 | 221.55 | 40 | 0 | 0 | 0 | 22.45 | 259.05 | 45 | 0 |
4-1430 | 0 | 0 | 13.53 | 250.89 | 42 | 0 | 0 | 0 | 22.16 | 368 | 41 | 0 |
4-1500 | 0 | 0 | 13.81 | 295.23 | 44 | 0 | 0 | 0 | 21.88 | 425.25 | 47 | 0 |
4-1530 | 0 | 0 | 15.7 | 441.1 | 45 | 0 | 0 | 0 | 20.34 | 408.15 | 46 | 0 |
4-1600 | 0 | 0 | 13.99 | 430.63 | 40 | 0 | 0 | 0 | 22.68 | 536.65 | 45 | 0 |
4-1630 | 105 | 0 | 138.64 | 697.47 | 45 | 0 | 105 | 0 | 53.45 | 461.14 | 53 | 0 |
4-1700 | 105 | 0 | 180.42 | 848.27 | 47 | 0 | 105 | 0 | 27.02 | 494.34 | 54 | 0 |
4-1730 | 105 | 0 | 62.14 | 878.77 | 53 | 0 | 105 | 0 | 27.16 | 1003.2 | 52 | 0 |
4-1800 | 105 | 0 | 125.7 | 710.72 | 48 | 0 | 105 | 0 | 25.36 | 427.58 | 43 | 0 |
4-1830 | 105 | 0 | 208.29 | 664.09 | 46 | 0 | 105 | 0 | 92.41 | 734.1 | 47 | 0 |
4-1900 | 105 | 0 | 133.17 | 493.23 | 49 | 0 | 105 | 0 | 63.86 | 338.05 | 43 | 0 |
4-1930 | 105 | 0 | 128.54 | 640.38 | 50 | 0 | 105 | 0 | 168.99 | 651.96 | 48 | 0 |
4-2000 | 105 | 0 | 90.8 | 507.54 | 47 | 0 | 105 | 0 | 65.67 | 556.58 | 50 | 0 |
4-2030 | 105 | 0 | 170.66 | 562.53 | 53 | 0 | 105 | 0 | 83.46 | 510.97 | 61 | 0 |
4-2100 | 0 | 0 | 12.62 | 569.4 | 49 | 0 | 0 | 0 | 21.21 | 376.63 | 45 | 0 |
4-2130 | 0 | 0 | 14 | 327.12 | 50 | 0 | 0 | 0 | 25.82 | 335.74 | 51 | 0 |
4-2200 | 0 | 0 | 15.35 | 314.91 | 51 | 0 | 0 | 0 | 22.69 | 404.91 | 57 | 0 |
4-2230 | 0 | 0 | 12.77 | 424.66 | 50 | 0 | 0 | 0 | 20.89 | 369.53 | 53 | 0 |
4-2300 | 0 | 0 | 13.5 | 510.89 | 75 | 0 | 0 | 0 | 21.67 | 228.67 | 52 | 0 |
4-2330 | 0 | 0 | 13.62 | 405.77 | 75 | 0 | 0 | 0 | 20.81 | 196.74 | 51 | 0 |
6-0000 | 0 | 0 | 13.43 | 214.68 | 58 | 0 | 0 | 0 | 20.37 | 247.61 | 43 | 0 |
6-0030 | 0 | 0 | 46.88 | 289.02 | 41 | 0 | 0 | 0 | 60.05 | 324.09 | 47 | 0 |
6-0100 | 0 | 0 | 159.84 | 464.4 | 39 | 1 | 0 | 0 | 133.52 | 496.28 | 39 | 0 |
6-0130 | 0 | 0 | 166.96 | 511.41 | 47 | 0 | 0 | 0 | 29.99 | 453.62 | 47 | 0 |
6-0200 | 0 | 0 | 1060.57 | 1398.41 | 44 | 1 | 0 | 0 | 29.12 | 389.04 | 40 | 0 |
6-0230 | 0 | 0 | 12.95 | 501.91 | 45 | 0 | 0 | 0 | 25.15 | 497.18 | 47 | 0 |
6-0300 | 0 | 0 | 40.01 | 488.95 | 46 | 0 | 0 | 0 | 29.78 | 550.22 | 44 | 0 |
6-0330 | 0 | 0 | 12.98 | 389.08 | 53 | 0 | 0 | 0 | 35.63 | 500.91 | 49 | 0 |
6-0400 | 0 | 0 | 31.08 | 421.51 | 48 | 0 | 0 | 0 | 32.45 | 353.1 | 51 | 0 |
6-0430 | 135 | 1 | 3600 | 3871.17 | 56 | 0 | 5 | 0 | 343.78 | 746.48 | 45 | 0 |
6-0500 | 3600 | 3600 | 5 | 0 | 423.29 | 781.24 | 45 | 0 | ||||
6-0530 | 135 | 0 | 2008.99 | 2334.06 | 52 | 0 | 135 | 0 | 1834.18 | 2210.58 | 51 | 0 |
6-0600 | 135 | 0 | 2162 | 2456.16 | 46 | 0 | 135 | 0 | 928.07 | 1209.29 | 47 | 0 |
6-0630 | 5 | 0 | 460.05 | 650.23 | 42 | 1 | 5 | 0 | 384.91 | 621.17 | 47 | 0 |
6-0700 | 5 | 0 | 2925.01 | 3250.83 | 42 | 0 | 5 | 0 | 337.99 | 631.76 | 33 | 2 |
6-0730 | 3600 | 3600 | 5 | 0 | 162.54 | 427.62 | 47 | 0 | ||||
6-0800 | 3600 | 3628.63 | 7 | 1 | 5 | 0 | 619.28 | 868.84 | 46 | 0 | ||
6-0830 | 3600 | 3600 | 5 | 0 | 403.02 | 658.92 | 51 | 1 | ||||
6-0900 | 130 | 0.96 | 3600 | 3950.85 | 50 | 0 | 5 | 0 | 182.53 | 503.37 | 48 | 0 |
6-0930 | 5 | 0 | 1745.26 | 2037.76 | 50 | 0 | 5 | 0 | 147.94 | 431.44 | 48 | 0 |
6-1000 | 5 | 0 | 517.31 | 904.9 | 47 | 0 | 5 | 0 | 172.34 | 487.23 | 44 | 0 |
6-1030 | 5 | 0 | 380.84 | 636.07 | 46 | 0 | 5 | 0 | 431.22 | 678.46 | 47 | 0 |
6-1100 | 5 | 0 | 1530.92 | 1762.32 | 56 | 0 | 5 | 0 | 2891.14 | 3193.12 | 49 | 0 |
6-1130 | 5 | 0 | 3295.1 | 3504.76 | 42 | 1 | 5 | 0 | 246.85 | 527.46 | 55 | 0 |
6-1200 | 0 | 0 | 1236.89 | 1604.46 | 65 | 0 | 0 | 0 | 78.82 | 438.08 | 66 | 0 |
6-1230 | 120 | 0 | 1810.36 | 2105.61 | 87 | 0 | 210 | 0.84 | 3600 | 3818.67 | 65 | 0 |
6-1300 | 120 | 0 | 3520.22 | 3854.08 | 56 | 0 | 120 | 0 | 184.92 | 477.78 | 69 | 0 |
6-1330 | 120 | 0 | 2996.05 | 3280.21 | 61 | 1 | 120 | 0 | 156.44 | 532.34 | 59 | 0 |
6-1400 | 3600 | 3600 | 120 | 0 | 630.71 | 872.7 | 69 | 0 | ||||
6-1430 | 120 | 0 | 1938.03 | 2284.72 | 61 | 0 | 120 | 0 | 1387.14 | 1600.86 | 54 | 1 |
6-1500 | 120 | 0 | 2044.11 | 2326.26 | 53 | 1 | 120 | 0 | 341.57 | 573.93 | 56 | 0 |
6-1530 | 10 | 0 | 1516.49 | 1997.61 | 95 | 0 | 10 | 0 | 1594.73 | 2090.87 | 60 | 0 |
6-1600 | 0 | 0 | 188.42 | 441.3 | 93 | 0 | 0 | 0 | 67.19 | 353.48 | 57 | 0 |
6-1630 | 0 | 0 | 15.09 | 285.28 | 53 | 0 | 0 | 0 | 23.73 | 299.96 | 53 | 0 |
6-1700 | 0 | 0 | 11.59 | 241.43 | 76 | 0 | 0 | 0 | 21.6 | 285.14 | 51 | 0 |
6-1730 | 0 | 0 | 12.23 | 200.72 | 63 | 0 | 0 | 0 | 23.41 | 245.59 | 51 | 0 |
INST | Val | Gap | AATM | PURE | SM-IT | NG | Val | Gap | AATM | HEUR | SM-IT | NG |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
6-1800 | 0 | 0 | 12.63 | 184.97 | 62 | 0 | 0 | 0 | 22.04 | 242.18 | 52 | 0 |
6-1830 | 105 | 0 | 30.04 | 253.18 | 53 | 0 | 105 | 0 | 28.19 | 217.35 | 53 | 0 |
6-1900 | 105 | 0 | 113.1 | 275.03 | 56 | 0 | 105 | 0 | 28.74 | 168.74 | 48 | 0 |
6-1930 | 105 | 0 | 24.55 | 162.25 | 52 | 0 | 105 | 0 | 30.25 | 336.35 | 113 | 0 |
6-2000 | 160 | 0 | 381 | 526.95 | 41 | 1 | 160 | 0 | 182.59 | 326.53 | 42 | 1 |
6-2030 | 105 | 0 | 24.6 | 206.29 | 54 | 0 | 105 | 0 | 30.56 | 214.95 | 56 | 0 |
6-2100 | 105 | 0 | 111.63 | 251.42 | 53 | 0 | 105 | 0 | 32.9 | 203.51 | 54 | 0 |
6-2130 | 105 | 0 | 131.22 | 289.15 | 54 | 0 | 105 | 0 | 95.31 | 242.89 | 54 | 0 |
6-2200 | 105 | 0 | 118.23 | 327.28 | 41 | 0 | 105 | 0 | 30.05 | 196.76 | 39 | 0 |
6-2230 | 105 | 0 | 115.56 | 318.96 | 62 | 0 | 105 | 0 | 29.13 | 198.5 | 42 | 0 |
6-2300 | 105 | 0 | 178.27 | 333.48 | 41 | 0 | 105 | 0 | 29.83 | 172.77 | 39 | 0 |
6-2330 | 105 | 0 | 219 | 405.63 | 39 | 0 | 105 | 0 | 31.44 | 185.82 | 38 | 0 |
7-0000 | 115 | 0 | 385.71 | 592.2 | 43 | 0 | 115 | 0 | 46.04 | 223 | 43 | 0 |
7-0030 | 115 | 0 | 286.23 | 485.91 | 38 | 1 | 115 | 0 | 56.98 | 203.08 | 37 | 1 |
7-0100 | 0 | 0 | 11.38 | 229.8 | 44 | 0 | 0 | 0 | 20.97 | 218.6 | 47 | 0 |
7-0130 | 0 | 0 | 13.79 | 266.62 | 42 | 0 | 0 | 0 | 27.9 | 266.67 | 61 | 0 |
7-0200 | 0 | 0 | 10.88 | 225.81 | 50 | 0 | 0 | 0 | 25.44 | 219.36 | 47 | 0 |
7-0230 | 0 | 0 | 11.06 | 208.87 | 48 | 0 | 0 | 0 | 41.09 | 216.97 | 45 | 0 |
7-0300 | 0 | 0 | 13.49 | 223.38 | 50 | 0 | 0 | 0 | 23.18 | 213.71 | 48 | 0 |
7-0330 | 0 | 0 | 12.28 | 216.05 | 51 | 0 | 0 | 0 | 25.18 | 231.28 | 52 | 0 |
7-0400 | 0 | 0 | 11.15 | 204.61 | 49 | 0 | 0 | 0 | 21.79 | 497.3 | 75 | 0 |
7-0430 | 0 | 0 | 13.11 | 244.61 | 63 | 0 | 0 | 0 | 42.55 | 632.78 | 104 | 0 |
7-0500 | 0 | 0 | 12.48 | 507.01 | 79 | 0 | 0 | 0 | 22.55 | 512.39 | 73 | 0 |
7-0530 | 0 | 0 | 11.08 | 551.25 | 88 | 0 | 0 | 0 | 22.81 | 376.79 | 65 | 0 |
7-0600 | 0 | 0 | 11.45 | 367.69 | 73 | 0 | 0 | 0 | 37.47 | 379.52 | 65 | 0 |
7-0630 | 0 | 0 | 13.1 | 390.19 | 74 | 0 | 0 | 0 | 21.98 | 271.49 | 52 | 0 |
7-0700 | 0 | 0 | 41.45 | 361.5 | 76 | 0 | 0 | 0 | 21.16 | 487.76 | 89 | 0 |
7-0730 | 0 | 0 | 11.52 | 400.81 | 87 | 0 | 0 | 0 | 21.69 | 443.46 | 98 | 0 |
7-0800 | 0 | 0 | 11.14 | 281 | 45 | 0 | 0 | 0 | 31.15 | 300.44 | 45 | 0 |
7-0830 | 0 | 0 | 12.09 | 263.23 | 52 | 0 | 0 | 0 | 40.3 | 296.65 | 43 | 0 |
7-0900 | 0 | 0 | 11.94 | 460.83 | 38 | 0 | 0 | 0 | 21.07 | 236.25 | 45 | 0 |
7-0930 | 0 | 0 | 11.19 | 389.66 | 38 | 0 | 0 | 0 | 21.12 | 204.68 | 46 | 0 |
7-1000 | 0 | 0 | 18.69 | 324.58 | 41 | 0 | 0 | 0 | 26.11 | 247.44 | 46 | 0 |
7-1030 | 0 | 0 | 13.76 | 275.46 | 48 | 0 | 0 | 0 | 44.57 | 264.75 | 45 | 0 |
7-1100 | 0 | 0 | 12.81 | 381.88 | 48 | 0 | 0 | 0 | 33.41 | 390.7 | 45 | 0 |
7-1130 | 0 | 0 | 25.76 | 562.52 | 48 | 0 | 0 | 0 | 28.36 | 417.74 | 55 | 0 |
7-1200 | 0 | 0 | 12.94 | 440.91 | 46 | 0 | 0 | 0 | 27.87 | 422.98 | 48 | 0 |
7-1230 | 0 | 0 | 235.79 | 741.71 | 79 | 1 | 0 | 0 | 71.63 | 335.89 | 70 | 0 |
7-1300 | 0 | 0 | 47.54 | 526.41 | 70 | 1 | 0 | 0 | 54.77 | 431.85 | 51 | 0 |
7-1330 | 0 | 0 | 12.08 | 359.71 | 57 | 0 | 0 | 0 | 59.25 | 403.07 | 54 | 0 |
7-1400 | 0 | 0 | 25 | 273.36 | 47 | 0 | 0 | 0 | 20.53 | 231.29 | 53 | 0 |
7-1430 | 85 | 0 | 761.41 | 951.07 | 47 | 0 | 85 | 0 | 325.23 | 589.92 | 58 | 0 |
7-1500 | 85 | 0 | 876.46 | 1168.02 | 53 | 0 | 85 | 0 | 73.89 | 328.41 | 48 | 0 |
7-1530 | 85 | 0 | 1756.19 | 2146.46 | 54 | 0 | 85 | 0 | 78.21 | 265.57 | 42 | 0 |
7-1600 | 120 | 0 | 3337.44 | 3521.73 | 56 | 0 | 120 | 0 | 101.14 | 297.24 | 55 | 0 |
7-1630 | 3600 | 3654.25 | 18 | 1 | 120 | 0 | 111.83 | 287.3 | 55 | 0 | ||
7-1700 | 120 | 0 | 1101.15 | 1299.68 | 54 | 0 | 120 | 0 | 104.34 | 298.16 | 55 | 0 |
7-1730 | 120 | 0 | 849.85 | 1026.91 | 53 | 0 | 120 | 0 | 94.22 | 286.12 | 55 | 0 |
7-1800 | 120 | 0 | 3241.61 | 3408.32 | 49 | 0 | 120 | 0 | 104.02 | 295.41 | 55 | 0 |
7-1830 | 120 | 0 | 1070.46 | 1273.82 | 56 | 0 | 120 | 0 | 129.86 | 332.25 | 58 | 0 |
7-1900 | 120 | 0 | 291.85 | 504.83 | 58 | 0 | 120 | 0 | 124.72 | 305.5 | 56 | 0 |
7-1930 | 120 | 0 | 209.13 | 400.24 | 53 | 0 | 120 | 0 | 160.53 | 329.07 | 54 | 0 |
7-2000 | 85 | 0 | 312.01 | 616.59 | 46 | 0 | 85 | 0 | 326.15 | 520.63 | 43 | 0 |
7-2030 | 0 | 0 | 11.96 | 240.67 | 50 | 0 | 0 | 0 | 24.21 | 205.72 | 42 | 0 |
7-2100 | 0 | 0 | 11.55 | 231.54 | 41 | 0 | 0 | 0 | 37.08 | 249.01 | 39 | 0 |
7-2130 | 0 | 0 | 12.03 | 253.55 | 64 | 0 | 0 | 0 | 21.26 | 311.16 | 43 | 0 |
7-2200 | 0 | 0 | 12.05 | 501.45 | 68 | 0 | 0 | 0 | 20.84 | 237.42 | 51 | 0 |
7-2230 | 3600 | 3612.98 | 4 | 2 | 3600 | 3615.93 | 4 | 3 | ||||
7-2300 | 0 | 0 | 160.01 | 555.43 | 42 | 0 | 0 | 0 | 20.63 | 596.2 | 60 | 0 |
7-2330 | 0 | 0 | 12.74 | 510.72 | 60 | 0 | 0 | 0 | 22.94 | 566.17 | 45 | 0 |
INST | Val | Gap | AATM | PURE | SM-IT | NG | Val | Gap | AATM | HEUR | SM-IT | NG |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
8-0000 | 0 | 0 | 12.54 | 431.94 | 59 | 0 | 0 | 0 | 20.02 | 393.66 | 56 | 0 |
8-0030 | 0 | 0 | 14.16 | 322.84 | 67 | 0 | 0 | 0 | 24.36 | 348.09 | 48 | 0 |
8-0100 | 0 | 0 | 13.6 | 228.38 | 57 | 0 | 0 | 0 | 24.3 | 375.52 | 50 | 0 |
8-0130 | 0 | 0 | 12.71 | 285.17 | 46 | 0 | 0 | 0 | 21.86 | 237.25 | 63 | 0 |
8-0200 | 0 | 0 | 12.31 | 211.97 | 54 | 0 | 0 | 0 | 21.08 | 339.5 | 49 | 0 |
8-0230 | 0 | 0 | 35.52 | 257.53 | 59 | 0 | 0 | 0 | 19.79 | 592.2 | 67 | 0 |
8-0300 | 0 | 0 | 12.35 | 227.02 | 63 | 0 | 0 | 0 | 18.75 | 253.74 | 60 | 0 |
8-0330 | 0 | 0 | 15.8 | 269.31 | 52 | 0 | 0 | 0 | 23.48 | 343.78 | 52 | 0 |
8-0400 | 0 | 0 | 12.58 | 264.62 | 56 | 0 | 0 | 0 | 23.59 | 454.84 | 51 | 0 |
8-0430 | 0 | 0 | 12.66 | 422.48 | 51 | 0 | 0 | 0 | 23.14 | 237.73 | 49 | 0 |
8-0500 | 0 | 0 | 15.33 | 202.67 | 49 | 0 | 0 | 0 | 21.08 | 373.8 | 63 | 0 |
8-0530 | 0 | 0 | 13.27 | 242.85 | 51 | 0 | 0 | 0 | 20.57 | 223.4 | 47 | 0 |
8-0600 | 0 | 0 | 11.28 | 201.32 | 51 | 0 | 0 | 0 | 19.78 | 264.12 | 47 | 0 |
8-0630 | 0 | 0 | 11.04 | 216.8 | 47 | 0 | 0 | 0 | 20.57 | 214.63 | 47 | 0 |
8-0700 | 0 | 0 | 144.53 | 369.6 | 44 | 1 | 0 | 0 | 19.74 | 236.69 | 48 | 0 |
8-0730 | 0 | 0 | 102.26 | 312.58 | 45 | 1 | 0 | 0 | 54.24 | 285.67 | 47 | 1 |
8-0800 | 0 | 0 | 25.16 | 312.84 | 78 | 1 | 0 | 0 | 73.82 | 340.47 | 43 | 1 |
8-0830 | 0 | 0 | 75.78 | 306.57 | 73 | 1 | 0 | 0 | 20.59 | 242.4 | 61 | 0 |
8-0900 | 0 | 0 | 178.06 | 416.84 | 58 | 1 | 0 | 0 | 19.87 | 321.84 | 60 | 0 |
8-0930 | 0 | 0 | 25.59 | 303.15 | 42 | 1 | 0 | 0 | 55.51 | 471.35 | 77 | 1 |
8-1000 | 0 | 0 | 54.47 | 432.24 | 55 | 1 | 0 | 0 | 57.19 | 242.46 | 38 | 1 |
8-1030 | 0 | 0 | 292.06 | 422.87 | 37 | 2 | 0 | 0 | 144.24 | 321.19 | 34 | 1 |
8-1100 | 0 | 0 | 26 | 325.5 | 64 | 1 | 0 | 0 | 21.97 | 224.77 | 58 | 0 |
8-1130 | 0 | 0 | 10.81 | 240.02 | 62 | 0 | 0 | 0 | 20.82 | 216.87 | 66 | 0 |
8-1200 | 0 | 0 | 12.25 | 253 | 72 | 0 | 0 | 0 | 22.63 | 190.99 | 50 | 0 |
8-1230 | 0 | 0 | 11.51 | 153.45 | 50 | 0 | 0 | 0 | 27.74 | 238.42 | 71 | 0 |
8-1300 | 0 | 0 | 11.88 | 175.78 | 54 | 0 | 0 | 0 | 25.62 | 214.49 | 55 | 0 |
8-1330 | 0 | 0 | 187.13 | 377.58 | 60 | 0 | 0 | 0 | 208.27 | 398.94 | 60 | 0 |
8-1400 | 0 | 0 | 12.66 | 190.35 | 45 | 0 | 0 | 0 | 88.46 | 283.25 | 62 | 0 |
8-1430 | 0 | 0 | 14.65 | 305.91 | 49 | 0 | 0 | 0 | 20.45 | 255.68 | 47 | 0 |
8-1500 | 0 | 0 | 13.82 | 312.95 | 46 | 0 | 0 | 0 | 20.1 | 184.6 | 42 | 0 |
8-1530 | 0 | 0 | 13.93 | 226.24 | 46 | 0 | 0 | 0 | 20.6 | 215.59 | 43 | 0 |
8-1600 | 0 | 0 | 11.83 | 233.74 | 52 | 0 | 0 | 0 | 20.74 | 246.94 | 51 | 0 |
8-1630 | 0 | 0 | 12.91 | 194.42 | 43 | 0 | 0 | 0 | 20.33 | 175.53 | 41 | 0 |
8-1700 | 0 | 0 | 12.84 | 348.89 | 45 | 0 | 0 | 0 | 21.51 | 171.66 | 42 | 0 |
8-1730 | 0 | 0 | 12.02 | 393.42 | 39 | 0 | 0 | 0 | 19.7 | 169.13 | 49 | 0 |
8-1800 | 0 | 0 | 12.04 | 283.66 | 41 | 0 | 0 | 0 | 19.91 | 191.79 | 60 | 0 |
8-1830 | 0 | 0 | 11.64 | 288.64 | 40 | 0 | 0 | 0 | 21.61 | 221.52 | 40 | 0 |
8-1900 | 0 | 0 | 12.67 | 208.95 | 43 | 0 | 0 | 0 | 18.88 | 205.89 | 43 | 0 |
8-1930 | 0 | 0 | 12.29 | 162.22 | 42 | 0 | 0 | 0 | 20.17 | 158.72 | 42 | 0 |
8-2000 | 0 | 0 | 39.09 | 203.99 | 41 | 0 | 0 | 0 | 19.09 | 223.35 | 44 | 0 |
8-2030 | 0 | 0 | 12.96 | 163.87 | 39 | 0 | 0 | 0 | 20.52 | 188.29 | 43 | 0 |
8-2100 | 0 | 0 | 12.09 | 274.03 | 39 | 0 | 0 | 0 | 17.73 | 170.16 | 48 | 0 |
8-2130 | 0 | 0 | 11.75 | 355.86 | 38 | 0 | 0 | 0 | 20.59 | 186.69 | 40 | 0 |
8-2200 | 0 | 0 | 14.32 | 261.25 | 39 | 0 | 0 | 0 | 21.18 | 199.56 | 45 | 0 |
8-2230 | 0 | 0 | 12.09 | 334.1 | 40 | 0 | 0 | 0 | 21.62 | 159.42 | 47 | 0 |
8-2300 | 0 | 0 | 13.24 | 244.29 | 48 | 0 | 0 | 0 | 21.81 | 299.87 | 38 | 0 |
8-2330 | 0 | 0 | 12.78 | 222.88 | 39 | 0 | 0 | 0 | 20.42 | 214.11 | 40 | 0 |
INST | Val | Gap | AATM | PURE | SM-IT | NG | Val | Gap | AATM | HEUR | SM-IT | NG |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
9-0000 | 0 | 0 | 13.65 | 309.57 | 43 | 0 | 0 | 0 | 19.48 | 323.47 | 43 | 0 |
9-0030 | 0 | 0 | 12.67 | 278.92 | 45 | 0 | 0 | 0 | 18.72 | 307.3 | 48 | 0 |
9-0100 | 0 | 0 | 13.94 | 296.77 | 46 | 0 | 0 | 0 | 23.55 | 436.77 | 43 | 0 |
9-0130 | 0 | 0 | 15.43 | 539.44 | 47 | 0 | 0 | 0 | 19.52 | 453.37 | 49 | 0 |
9-0200 | 0 | 0 | 12.18 | 435.51 | 47 | 0 | 0 | 0 | 20.14 | 393.52 | 47 | 0 |
9-0230 | 0 | 0 | 14.17 | 388.23 | 48 | 0 | 0 | 0 | 20.4 | 271.1 | 45 | 0 |
9-0300 | 0 | 0 | 14.27 | 393.66 | 46 | 0 | 0 | 0 | 20.02 | 261.25 | 46 | 0 |
9-0330 | 0 | 0 | 17.09 | 545.31 | 50 | 0 | 0 | 0 | 19.41 | 426.01 | 47 | 0 |
9-0400 | 0 | 0 | 15.21 | 441.35 | 47 | 0 | 0 | 0 | 19.81 | 319.63 | 47 | 0 |
9-0430 | 0 | 0 | 21.13 | 457.47 | 48 | 0 | 0 | 0 | 21.39 | 314.79 | 47 | 0 |
9-0500 | 0 | 0 | 12.94 | 248.84 | 47 | 0 | 0 | 0 | 25.42 | 358.44 | 47 | 0 |
9-0530 | 0 | 0 | 15.13 | 403.15 | 47 | 0 | 0 | 0 | 21.33 | 467.63 | 49 | 0 |
9-0600 | 0 | 0 | 16.19 | 272.61 | 47 | 0 | 0 | 0 | 21.02 | 257.6 | 48 | 0 |
9-0630 | 0 | 0 | 18.99 | 311.4 | 46 | 0 | 0 | 0 | 20.12 | 339.96 | 46 | 0 |
9-0700 | 0 | 0 | 14.5 | 325.94 | 46 | 0 | 0 | 0 | 20.79 | 361.05 | 44 | 0 |
9-0730 | 0 | 0 | 16.44 | 381.39 | 55 | 0 | 0 | 0 | 19.82 | 334.4 | 45 | 0 |
9-0800 | 0 | 0 | 14.35 | 319.69 | 43 | 0 | 0 | 0 | 19.31 | 421.98 | 45 | 0 |
9-0830 | 0 | 0 | 14.91 | 367.42 | 43 | 0 | 0 | 0 | 97.13 | 347.54 | 34 | 2 |
9-0900 | 0 | 0 | 14.9 | 365.14 | 43 | 0 | 0 | 0 | 60.66 | 343.82 | 39 | 1 |
9-0930 | 0 | 0 | 17.45 | 255.33 | 45 | 0 | 0 | 0 | 60.27 | 268.48 | 40 | 1 |
Rights and permissions
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 Commons 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://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Hoppmann-Baum, K., Hennings, F., Lenz, R. et al. Optimal Operation of Transient Gas Transport Networks. Optim Eng 22, 735–781 (2021). https://doi.org/10.1007/s11081-020-09584-x
Received:
Revised:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s11081-020-09584-x