Optimizing transient gas network control for challenging real-world instances using MIP-based heuristics

Optimizing the transient control of gas networks is a highly challenging task. The corresponding model incorporates the combinatorial complexity of determining the settings for the many active elements as well as the non-linear and non-convex nature of the physical and technical principles of gas transport. In this paper, we present the latest improvements of our ongoing work to solve this problem for real-world, large-scale problem instances: By adjusting our mixed-integer non-linear programming model regarding the gas compression capabilities in the network, we reflect the technical limits of the underlying units more accurately while maintaining a similar overall model size. In addition, we introduce a new algorithmic approach that is based on splitting the complexity of the problem by first finding assignments for discrete variables and then determining the continuous variables as locally optimal solution of the corresponding non-linear program. For the first task, we design multiple different heuristics based on concepts for general time-expanded optimization problems that find solutions by solving a sequence of sub-problems defined on reduced time horizons. To demonstrate the competitiveness of our approach, we test our algorithm on particularly challenging historic demand scenarios. The results show that high-quality solutions are obtained reliably within short solving times, making the algorithm well-suited to be applied at the core of time-critical industrial applications.


Introduction
The natural gas supply of Europe and especially Germany is subject to a disruptive transformation.With the possibility of cutting the supply from Russia and the recently solidified plans to build LNG terminals in Germany as published by the government [Federal Ministry for Economic Affairs and Climate Action, 2022], new flow patterns are likely to emerge in the corresponding transport networks.Moreover, the usage of fossil-based natural gas will soon be terminated in favor of more sustainable alternatives as part of the energy transformation towards massive reductions in CO 2 emissions.In particular, the natural gas networks may be used to store synthesized methane generated via power-to-gas from excess renewable energy or could (partially) be repurposed for hydrogen transport, see, for example, the scenario framework for the German Gas Network Development Plan as published by the German gas network operators [FNB Gas -Association of German transmission system operators, 2021].Especially when applying power-togas strategies, the intraday supply patterns are likely to become more volatile, complicating the control of the networks even further.
Mathematical optimization can assist the network operators in preparing for these new challenges and ensuring a secure and efficient network operation.However, solving the corresponding optimization problems is hard.On the one hand, the physics of gas transport in pipelines describing the interplay of gas pressure and flow is captured by non-convex equations.On the other hand, the network contains many active elements to control the gas flow, like valves to dynamically change the network topology and compressors to increase the gas pressure.Apart from the non-convex description of the maximum power available to compressors, these control options also introduce an immense combinatorial complexity to the problem.Since the problem is both practically relevant and theoretically challenging, it has been intensively studied in the past, and we will give an overview split into two categories in the following.
First, there is the stationary variant of the problem.It defines an equilibrium state of the network by assuming a single infinitely long time step and can, for example, be used for planning purposes.For a profound introduction to stationary gas network optimization, we recommend the book of [Koch et al., 2015] or the paper of [Pfetsch et al., 2015] for a summary of the corresponding findings.They model the problem as a mixed-integer non-linear program (MINLP) and deal with the two sources of complexity separately by first finding solutions for the discrete decision variables and afterwards completing these to highly precise solutions by solving a corresponding non-linear programming (NLP) model.A more recent study on the stationary problem is presented in [González Rueda et al., 2019].They use a model based on simplified structures representing the network areas containing the compressors.Algorithmically, the authors mainly focus on solving an NLP-model variant of the problem with a two-step sequential linear programming (SLP) approach, which is inspired by an early attempt at solving the problem reported in [Pratt and Wilson, 1984].However, by replacing the first step with sequential mixed-integer programming, they can convert their algorithm into a heuristic for MINLP problems.Test on real-world instances, which only have a limited number of binary variables, reveal a performance similar to state-of-the-art NLP and MINLP solvers.
Second, the transient variant of the problem aims at finding a short-term control for the network over discrete future time steps.An early publication on the topic was the thesis of [Moritz, 2007].
Here, the problem was formulated as a pure mixed-integer program (MIP) by approximating the nonlinear pipe equations and compression power constraints using piece-wise linear functions.To aid the corresponding solver, the author proposed a special branching scheme for the piece-wise linear functions and used a simulated annealing heuristic, which was also described in [Mahlke et al., 2007].Shortly after, [Domschke et al., 2011] presented their approach for a similar problem formulation, which again used piece-wise linear functions to approximate the non-linear parts of the network.They iteratively improved the linearization by first solving the corresponding MIP model, then using the binary solution values to form a non-linear programming model using the original not approximated equations, and afterwards refining the linearization by adding mesh-points corresponding to the found solutions.The procedure is stopped when a predefined approximation threshold is satisfied.In their computational results, they report problems in solving the single MIP models to optimality even after multiple hours.In a more recent study, [Hahn et al., 2017] directly solve an MINLP formulation of the problem.Although the problem is non-convex, they use a convex solver and apply a customized branching rule for the problem.In their experiments, they design scenarios for small test networks of less than 10 nodes that include the reversal of the gas flow direction on single pipelines, which was not possible in the networks used in [Moritz, 2007] and [Domschke et al., 2011].Despite the size of the networks, the different variants of their algorithm need, even in the best cases, hours to finish.Another contribution was made by [Burlacu et al., 2019].By introducing a new discretization for the differential equations describing the gas flow in pipes, they were able to reformulate the transient problem as a stationary one on a time-expanded graph and use a specialized algorithm designed for the stationary case to solve the problem.It is based on iteratively solving mixed-integer programming relaxations of the problem, which they again obtained by applying a piece-wise linearization.By refining the linearization in each iteration, their approach converges to the globally optimal solution.To speed up the solving process, the intermediate MIP solutions were completed to valid primal solutions of the overall problem using an NLP solver.For the one evaluated instance, the algorithm finds a solution with a gap of less than 10% in less than 60 seconds and less than 8% in half an hour.However, after this point, even single iterations of the algorithm reach run times of over an hour and may not even improve the gap.
Finally, we mention the work of [Hoppmann-Baum et al., 2021], who presented a decision support system for network operators to ensure a secure and efficient network operation.The system is split into two separate stages: First, an MINLP problem formulation covering the complete network is solved, which approximates the capabilities of the technical elements located at large pipeline intersections, the network stations.In particular, each of these stations is replaced with a simplified graph model using artificial arcs.Non-technical control measures available to the network operators are realized as slack variables to adjust the future supplies and demands in terms of pressure and inflow.Their usage as well as changes in the technical control are penalized in the objective functions of a tri-level optimization model.To derive a feasible solution, a series of linearized MIP models are solved, which use two types of primal heuristics to determine good initial solutions.Afterwards, a linear programming (LP) model is solved to smooth the solution, followed by an SLP approach that transforms the solution into one respecting the original non-linear constraints of the MINLP problem.In a computational study on subsequent real-world instances spanning multiple days, the algorithm's run time and the differences in the solutions are evaluated.In the second stage of the decision support system, the actual control measures recommended to the network operators are determined by solving highly detailed models for each network station individually, as described in [Hennings et al., 2021a].Note that since the overall system runs in regular intervals to provide up-to-date control recommendations, it must reliably satisfy strict time limit restrictions.
In this paper, we present several improvements to the simplified network-wide model of the first stage of the decision support system described in [Hoppmann-Baum et al., 2021]: We suggest an alternative model for the artificial network station arcs representing the gas compressors, which is closer to the compression capabilities of the highly detailed network station model used in the second stage without increasing the model's size.Furthermore, we formulate a new algorithm to solve the MINLP problem of the first stage, which finds better solutions faster than the original approach.In particular, we first determine sets of values for all the binary variables of the model, which we call binary assignments, and afterwards complete these to full solutions of the MINLP by solving the remaining non-linear programming (NLP) model.The binary assignments are obtained by combining a sequential mixed-integer programming (SMIP) routine with heuristics concepts for time-expanded problems: The rolling horizon and the aggregate horizon heuristic.We additionally apply a dynamic branch-and-bound node limit to the single model solves to further reduce the maximum run times.Finally, we thoroughly test the proposed algorithm's run time and its solution quality against valid lower bounds.As a test set, we identify especially hard-to-solve realworld instances of the problem that feature many supply and demand changes.Thereby we verify that our algorithm is suitable for time-critical industry applications and more volatile future supply patterns.
Our proposed approach focuses on finding good primal solutions fast and does not provide a valid lower bound on the value of the problem's optimal solution.This enables us to tackle larger and more challenging problem instances than those found in the literature while respecting strict time limits.We follow the idea of dealing with the combinatorial and non-linear aspects of the problem separately and also use sequential mixed-integer programming for finding good binary assignments.However, since we do not aim at finding lower bounds, we do not need to use an iteratively refined piece-wise linear formulation but only an updatable linearization, similar to what [González Rueda et al., 2019] propose for the first level of their algorithm.As a consequence, the general size of the MIP problems is smaller and does not increase during the solving process, which is one of the reasons for the reported slow lower bound convergence of the transient approaches listed above.
We note that the heuristic concepts we use for time-expanded problems are well known in the literature.The rolling horizon approach was, for example, used for the optimization of product manufacturing in [Tiacci and Saetta, 2012], distributed power generation in [Capitanescu, 2017], or as one of the primal heuristics in [Hoppmann-Baum et al., 2021].A theoretical analysis of rolling horizon approaches, as well as a more extensive list of applications, is given in [Glomb et al., 2022].Similarly, heuristics using an aggregation of a time horizon were often described, as in the example of scheduling the iron ore production in [Newman and Kuchta, 2007] or the extension of freight transport networks in [Boland et al., 2013].
The rest of the paper is structured as follows: Section 2 introduces the MINLP problem formulation, where we particularly highlight the differences with respect to the model used in [Hoppmann-Baum et al., 2021].In the following Section 3, we present our heuristic algorithmic approach for solving this problem.The computational results are the content of Section 4, where we describe the creation of the instance sets, compare a base version of our algorithm against a global MINLP solver in terms of solution quality on smaller instances, and then compare the actually proposed algorithm against the base algorithm version on a more challenging instance set.We close with concluding remarks in Section 5.

Model formulation
Our problem is based on the gas flow model used for the first stage of the decision support system presented in [Hoppmann-Baum et al., 2021].As in there, we assume that the gas network contains a set I of network stations.Those are sub-networks that include the majority of all actively controllable network elements, among them all the network's compressor stations, and are mainly located at large pipeline intersections.For the model, the inner network topology of each station is replaced by a set of artificial arcs representing an approximation of the station's gas routing, compression, and regulation capabilities.This artificial model was created manually by industry experts who work for the network operator and have an excellent knowledge of the individual stations and their properties.
In the following, we summarize the differences between our problem formulation and the one of Hoppmann-Baum et al..Note that we refer to their model as the base model for the remainder of the paper.
• As the most crucial change, we introduce in Section 2.5.4 a new model for the artificial compressor arcs in the network stations.It is based on the usage of configurations for each artificial arc, which we construct based on the properties of the compressor stations of the original network station topology.We show that the model is much closer to the highly detailed compressor station model used in [Hennings et al., 2021a] while maintaining a similar model size with respect to the number of variables and constraints.
• The second major difference is the introduction of artificial valves at the boundaries of the network stations in Section 2.7.With these, it is possible to decouple the pressure inside and outside the station at these points.Note that this was initially suggested in [Hoppmann-Baum, 2022].
• In this paper, we do not use a tri-level model to reflect strict objective function priorities for using certain gas network control measures.Instead, we use a single objective function representing the weighted sum of different terms to penalize them, see Section 2.10.This objective also includes smoothing terms, which minimize pressure and inflow changes at the network station boundaries over time.As a final change in the objective, we include a newly created incentive to reduce the usage of compressor units based on the new artificial compressor configurations to avoid unnecessary energy consumption.
• We use the regulator model of [Hennings et al., 2021a], which is slightly different from the base regulator model.
• In the artificial arc models introduced in Section 2.5, we do not include combined arcs, which have gas compression as well as regulation capabilities, as they do not appear in our real-world-based test instances, see Table 2 in Section 4.1.For the same reason, we did not include artificial bi-directional compressor arcs or bi-directional combined arcs.
• Finally, we also do not introduce flow direction conditions as they are not used in the network station description of our test instances.
When formulating our model below, we often use indicator constraints, which have two parts: The first one is a logical condition in the form of a binary variable y attaining a value b ∈ {0, 1} and the second one an implied linear constraint a T x ≤ a 0 of variables x = {x 1 , . . ., x n } and coefficients a ∈ R n , a 0 ∈ R, which only holds if the logical condition is true.We denote such an indicator constraint by Note that we also use as an abbreviation of the two indicator constraints Indicator constraints can be reformulated as linear constraints via the big-M approach in case all the variables x = {x 1 , . . ., x n } are bounded, see [Bonami et al., 2015] for a general introduction to the topic.
A list of all the variables used in our model, together with their domains, meanings, and units, is given in Table A.1 in the Appendix.

Gas flow in networks
We model a gas network as a directed graph G = (V, A) with V being the set of nodes or vertices and A = A pi ∪ A va ∪ A rg ∪ A ar being the set of arcs representing the single element sets present in the network: Pipes A pi , valves A va , regulators A rg , and artificial network station arcs A ar .The latter is again composed of the different types of artificial arcs found in the network stations: Shortcuts A arSc , regulating arcs A arRg , compressor arcs A arCp , and bi-regulating arcs A arBiRg .For station i ∈ I, we denote by we denote its set of artificial arcs and by V fn i ⊆ V the nodes forming the station's boundary, which we call fence nodes.Formally, the fence nodes are those nodes having incident arcs from within the station as well as from outside of the station, i.e., ∀v ∈ V fn i : (∃a in ∈ A ar i : Note that, in theory, a node can be the fence node of different stations.However, this does not happen in practice.Furthermore, there can be artificial non-fence nodes inside the station, which are needed for the manual station modeling.Finally, we assume for each station i ∈ I that the graph induced by the station's arcs A ar i is connected.Furthermore, we consider a discrete time horizon T 0 := {0, . . ., t max }, where T := T 0 \ {0} denotes the set of all future time steps.We denote by τ (t) the time difference in seconds from timestep t ∈ T 0 to the initial state time step 0.
The gas flow on each arc a = ( , r) from node to node r at some time t ∈ T 0 is characterized by the pressure p and the mass flow q at each end of the arc, i.e., the incoming pressure p ,a,t , the outgoing pressure p r,a,t , the incoming mass flow q ,a,t , and the outgoing mass flow q r,a,t .Here, negative flow values represent flow against the arc's orientation, while the pressure values are always positive.Since we assume that only pipes have a nonzero volume, the inflow and outflow of all non-pipe arcs is identical.Thus, we define the mass flow of an arc a ∈ A \ A pi as q a,t := q ,a,t = q r,a,t .
The quantities of the arcs are connected at the nodes of the network.For all arcs incident to a node, the corresponding pressure values at that end node are equal.Hence, we can reference this value as pressure p v,t of node v ∈ V for time t ∈ T 0 and it holds: For the mass flow holds the conservation of mass, a Kirchhoff-type law which demands that the flow into and out of a node v ∈ V must be balanced for all time steps t ∈ T , i.e., a∈A:a=(v,r) (1) This equation is also known as flow balance equation.Note that the inflow d v,t denotes gas injection into node v at time t from outside of the network for positive values and gas withdrawal for negative values.Furthermore, the inflow value at time t represents the average gas inflow rate into the node in the continuous time interval between t-1 and t.We partition the nodes into three sets: The set of entries , and the set of inner nodes V 0 with d v,t = 0 for v ∈ V 0 for all t ∈ T .We call the union of all entries and exits the boundary nodes V b of the network.
For each pressure and flow variable, we are given lower and upper bounds for each t ∈ T , which are denoted by bars below and above the corresponding variable.
Apart from the network topology, we are also given the initial state of the network, which specifies for t = 0 the node pressure p v,0 for nodes v ∈ V and the arc flow values q l,a,0 and q r,a,0 for arcs a = (l, r) ∈ A. Furthermore, the initial states of the active network elements outside the network stations, i.e., valves and regulators, are given.

Pipes
For pipes, we use the friction-dominated model given as (ISO3) in [Domschke et al., 2021].It is based on the Euler equations modeling the one-dimensional gas flow in straight cylindric pipes and can be derived by assuming a time-constant temperature (isothermal conditions), gas mixture, and compressibility factor, as well as clearly subsonic gas velocities and a negligible inertia term, which both are usually given in real-world conditions, see [Osiadacz, 1996;Domschke et al., 2021] and [Hennings, 2021].The usage of the implicit box scheme of [Domschke et al., 2011;Kolb et al., 2010] as a discretization yields the model for a pipe a = ( , r) ∈ A pi and two adjacent time points t 1 ∈ T and t 0 := t 1 − 1 given as The first equation (2a) is called the Continuity equation, while the second equation ( 2b) is known as the Momentum equation and contains the non-linear friction term.In the equations, the specific gas constant is given by R s , the gravitational acceleration by g, the gas temperature of the pipe by T a , its compressibility factor by z a , its length by L a , its diameter by D a , its area by A a defined as A a = D 2 a π/4, its Darcy friction coefficient by λ a , and finally its slope by s a defined as s a = (h r − h )/L a using the height or altitude h of the pipes' end nodes.The gas temperature depends on the individual pipe, as we determine it based on the given gas temperatures T v,0 in the initial state for each node v as In a similar fashion, we determine the pipes' compressibility factor using the formula of Pápay [Pápay, 1968;Takacs, 1989] based on the initial pressure and temperature of each node as z a = 0.5(z(p ,0 , T ,0 ) + z(p r,0 , T r,0 )).
Finally, the Darcy friction factor is determined using the formula of Nikuradse [Nikuradse, 1950], which is an explicit approximation depending only on the pipe's diameter and its integral roughness k a .

Valves
Valves are active elements that are able to dynamically connect or disconnect two network nodes.Their state is represented by a binary mode variable m op a,t , which specifies if the valve a is open (m op a,t = 1) at time t, connecting both end nodes and forcing their pressures variables to attain identical values, or closed (m op a,t = 0), disconnecting the nodes, decoupling the respective pressure values, and prohibiting gas flow exchange.Thus, the corresponding model for a valve a = ( , r) ∈ A va reads

Regulators
Regulators are active elements that can be seen as an extension of valves.In addition to changing the network's connectivity, they can reduce the pressure of the gas in the direction of the flow.However, this is only possible along the arc's orientation.We use the regulator model of [Hennings et al., 2021a] Note that in the indicator condition of Equation (4c), the term m ac a + m op a acts as a binary variable, as Equation (4a) forces the choice of exactly one of the three modes.Furthermore, we assume an internal check valve in the element, which prevents flow against the orientation even in the case of an open regulator.For a more detailed description of the internal regulator behavior, we refer to [Hennings et al., 2021b].

Artificial station arcs
As stated above, the original network topology inside the stations is replaced by artificial arcs of different types, which represent an approximation of the station's gas transport capabilities.Each arc can be active or inactive, which we model by a binary variable where x arc a,t = 1 represents the active state.An inactive arc behaves like a closed valve, i.e., the flow is zero and the end nodes' pressures are decoupled, while the active behavior is specific to each artificial arc type, which we discuss in the following.

Shortcuts
Shortcuts connect its end nodes if active.Hence, their model is equivalent to that of a valve and can be stated for a shortcut a = ( , r) ∈ A arSc and time t ∈ T as: x arc a,t = 0 =⇒ q a,t = 0. (5b)

Regulating arcs
Artificial regulating arcs capture the general behavior of regulators.Hence, they only allow flow along their orientation and offer the possibility of pressure reduction if they are active.The corresponding model for a regulating arc a = ( , r) ∈ A arRg and time t ∈ T can be stated as x arc a,t = 0 =⇒ q a,t ≤ 0 (6b) Note that the Equations (4) do not apply here.

Bi-regulating arcs
To represent the behavior of a pair of anti-parallel regulators in the original network topology, we use a bi-regulating arc.We can model its behavior by choosing an orientation of the element in active mode and applying the corresponding constraints of a regulating arc for each direction separately.The model for a bi-regulating arc a = ( , r) ∈ A arBiRg and time t ∈ T is given as q a,t = q fwd a,t − q bwd a,t x fwd a,t = 0 =⇒ q fwd a,t ≤ 0 (7d) q fwd a,t ≥ 0 (7e) x bwd a,t = 0 =⇒ q bwd a,t ≤ 0 (7g) q bwd a,t ≥ 0 (7h) x fwd a,t , x bwd a,t ∈ {0, 1}.
Here, x fwd a,t represents the choice for the forward orientation for arc a and time t, in which flow and regulation happen along the arc's orientation, while x bwd a,t represents the choice of the backward orientation, in which the orientation is flipped.In the same fashion, q fwd a,t represents the forward flow while q bwd a,t represents the flow in the backward direction.

Compressor arcs
Finally, there are the compressor arcs, which are the only elements able to increase the gas pressure in the direction of the flow.They approximate the capabilities of compressor station arcs in the original network station topology.For this reason, we shortly introduce these in the following.Each compressor station, which we assume to be modeled as an arc a = ( , r) as well, represents a combination of single compressor units, which again are a combination of a compressor machine and a corresponding drive that provides the necessary power.The compressor machine of a compressor unit u can operate in its feasible operating range, defined as a polytope in the space (Q, ψ) of the volumetric flow Q and the compression ratio ψ given as Here, ρ denotes the density of the gas, ρ the incoming gas density, while z a and T a are defined in the same way as they were for pipes in Section 2.2.Equation (8a) is called the equation of state for real gases [Osiadacz, 1996;Fügenschuh et al., 2015].Furthermore, the compression is restricted by the maximal amount of power Pu the corresponding drive can provide.The power needed for compressing a mass flow value of q from p to p r , i.e., with ψ = p r /p , is given as Here, η ad is the adiabatic efficiency and is assumed to be a constant given per compressor unit.Further, κ is the isentropic exponent, which we approximate by the constant value of 1.296, see, for example, [Fügenschuh et al., 2015] for more details.The single units can then be combined to predefined configurations C a of the compressor station a, defined as serial connections of parallel unit arrangements.For a more in-depth description, we refer to [Fügenschuh et al., 2015;Walther and Hiller, 2017;Hennings et al., 2021a].
In the base model of [Hoppmann-Baum et al., 2021], the artificial pressure increasing arcs were modeled as a single compressor representing the combination of single units, which could dynamically be assigned to pressure-increasing arcs using binary assignment variables.As a consequence of this dynamic combination of different units, some of the constraints describing the resulting compression capabilities were approximated quite generously to keep the overall model linear.
To avoid this problem, we develop a novel formulation, which is closer to the original compressor station model while the model size with respect to the number of variables and constraints is not increased.Like for the compressor stations, we create for each artificial compressor arc a = ( , r) ∈ A arCp a set of configurations C ar a to choose from.These form different compression settings based on how the compressor units can be composed in the underlying compressor stations.The feasible operating range of a chosen configuration c ∈ C ar a of an artificial compressor arc a is restricted in terms of three parameters: The maximum compression ratio ψ c , the maximum volumetric flow Q c , and the maximum compression power P c .Furthermore, each configuration has an associated number u c of used compressor units.From the parameters restricting the feasible operating range, we further find parameters α p c , α pr c , and α q c , which represent a linearization of the compression power of the form and are determined by a linear regression based on triples (p s , p s r , q s ) s∈N sampled from within the intersection of the corresponding variable bound intervals over time, which satisfy the constraints Note that usually, α p c < 0, α pr c > 0, and α q c > 0 hold.The final component of the configuration-based model is the set of artificial compressor configuration conflicts Z.Each conflict z is defined as tuple z = (c 1 , c 2 ) of configurations of different artificial compressors from the same network station that cannot be used simultaneously.Conflicts are needed to model, for example, compressor units that can be used in different compressor stations, but not at the same time.As a consequence, configurations c 1 and c 2 of two corresponding artificial compressors cannot be used simultaneously if both can use the same unit.Since conflicts are always defined for configurations of artificial compressors in the same network station, we can partition the conflicts per network station as Z = ˙ i∈I Z i .Given this description, we state the configuration-based model for an artificial compressor a = ( , r) ∈ A arCp and time t ∈ T as x arc a,t = 0 =⇒ q a,t ≤ 0 (10c) where we additionally demand for each configuration c ∈ C ar a of an artificial compressor a = ( , r) Assume a artificial compressor a ac with a single associated compressor station a cs .
Let the configurations of a cs be defined as C acs := .
Then the corresponding artificial configurations for a ac are C ar aac := . 1 , 1 , {u 1 , u 2 }),  ( 2 , 1 , {u 1 , u 3 }),  ( 2 , 1 , {u 2 , u 3 }) Figure 1: Example of the construction of the configurations for an artificial compressor a ac .The associated compressor station a cs has four different configurations using three different compressor units: Two configurations using only a single unit, either unit u 1 or unit u 2 , and two configurations in which each of the two units is used in parallel to unit u 3 .There are three resulting artificial configurations for a ac , which all allow only one serial stage: One configuration with only one allowed parallel unit, in which u 1 or u 2 can be chosen, and two configurations allowing for two parallel units, one with the units set {u 1 , u 3 } and one with the set {u 2 , u 3 }.
While we can combine the first two original configurations into one artificial configuration, this is not possible for the parallel configurations.The corresponding unit set would have to be {u 1 , u 2 , u 3 }, which is not feasible as there is no original configuration combining u 1 and u 2 in parallel. that Here, the binary variable x cfg c,t represents the decision to use (value=1) or not to use (value=0) the corresponding configuration c ∈ C ar a at time t.In addition, we need the following constraint for each conflict z = (c 1 , c 2 ) ∈ Z and time t ∈ T : Note that the artificial compressor arcs are directed and therefore can only have positive flow as well as compress gas from the end node towards r, which is ensured via Equations ( 10d) and (10b).
Artificial configuration set construction Since the artificial compressor configurations are not given as part of the input, we create them for each artificial compressor a ∈ A arCp based on the original compressor stations represented by it.Here, we associate each artificial configuration c with a tuple (u p c , u s c , U c ), where u p c is the number of parallel compressor units to use, u s c is the number of serial units to use, and U c is the set of usable units.We manually create these tuples based on the configurations of the original compressor stations as well as their topological position in the network.If multiple original compressor compositions use similar units in the same layout, we represent them by the same artificial configuration if possible.An example of the construction is given in Figure 1.
Given (u p c , u s c , U c ) for each artificial configuration c, we now derive the corresponding parameters for the number of simultaneously used compressor units u c , the maximum compression ratio ψ c , the maximum volumetric flow Q c , and the maximum power P c .First, we present a procedure for a configuration c for the special case of either u p c or u s c being equal to 1, i.e., for the case of a pure serial or pure parallel compression.For each compressor unit u ∈ U c , we are given the maximum power Pu its associated drive can provide.Furthermore, we can determine its maximum compression ratio ψu and maximum volumetric flow Qu by iterating the extreme points of its feasible operating range polytope.Let then ( ψn ) |Uc| n=1 be the Figure 2: Relative error in compression power.We compare the power P sum needed for a serial compression using two separate compression processes with ratios ψ 1 and ψ 2 to the power P one needed by our approximation of using a single compression process with ratio ψ 1 • ψ 2 .We plotted the relative error (P one − P sum )/P sum for compression ratios ψ 1 and ψ 2 up to 2.0.
Note that for the serial compression, summing up the power values is not quite accurate and slightly underestimates the necessary compression power when comparing the original two serial compression processes with one combined compression process using the product of the single compression ratios, see Figure 2. As a consequence, P c is potentially too small to allow for all compression processes possible in the original network.However, in the data used for our computational experiments in Section 4, the maximum serial compression number u s c for a configuration c is 2. Furthermore, all involved units have compressor ratios smaller than 2.0.Thus, the approximation error is guaranteed to be smaller than 9%.
If a configuration features multiple serial and multiple parallel stages, the above-presented formulas for creating the configuration parameters cannot be applied directly.Instead, we create for each such arrangement and each concrete assignment of a compressor unit to a position in the arrangement a separate artificial configuration for the corresponding artificial compressor arc.The total number of used compressor units u c is then equal to the number of active units in each specific arrangement, and the maximum power P c is determined as the sum of maximum power values of these units.To determine the parameters ψ c and Q c , we recursively apply the formulas used in Equations ( 13b) and (13c) for each of the concrete serial and parallel sub-compositions of units.The approach of creating one configuration for each possible unit arrangement can, in theory, result in a lot of artificial configurations.However, in our experience, this is not the case in practice, as such configurations do not occur very often.

Model comparison
The main benefit of the new artificial compressor model in Equations ( 10), (11), and ( 12) compared to the base model is that each compressor unit contributes either in serial or in parallel to the compression capabilities and does not use the benefits both at the same time.
Furthermore, we improved the description of the maximum compression ratio ψ and the maximum volumetric flow Q.For both, we do not the need to replace the incoming pressure of the arc p by a constant value, which was done explicitly in Equation ( 40) in [Hoppmann-Baum et al., 2021] and implicitly by assuming a maximum mass flow bound per compressor unit, which can only be derived from the inherent volumetric flow bound by assuming a constant incoming pressure.Furthermore, we can correctly employ the product of compressor ratios in the computation of the maximum usable compression ratio of serial compression processes in Equation (13b) instead of the linear approximation formula used in ( 35) in [Hoppmann-Baum et al., 2021].
As an example, we show the model improvements in comparison to [Hoppmann-Baum et al., 2021] for a single artificial compressor arc in Figure 3, where we compare the feasible region polytopes produced by the base model, our new formulation, and the actual feasible region of the configuration polytopes of the compressor station in the original station topology.To have a clean comparison, we choose an artificial compressor arc representing exactly one compressor station.This compressor station features three compressor units, from which two units can be combined in parallel, but has no serial unit combination capabilities.Note that we excluded the power bound from the feasible region polytopes since it is modeled for both models in the same way, i.e., as a linear approximation using a linear regression approach.
We clearly see the increased accuracy of the new model already for an artificial compressor arc, which is only able to combine two compressor units.When combining more compressor units, the differences would further increase.

Station simple states
The stations' artificial arcs do not operate independently from each other.Instead, each station i ∈ I is always in exactly one simple state, where the set of possible simple states of the station is given as S i , and the union of simple states over all stations is given as S = i∈I S i .Each simple state s ∈ S i restricts the usage of the artificial arcs by defining a set of arcs A on s ⊆ A

Station flow directions and fence node valves
Apart from the simple state, each station is also assigned exactly one flow direction f at each point in time t ∈ T 0 , captured by the binary variable x fd f,t , where the value 1 again stands for an active flow direction.The set of possible flow directions for station i is given as F i .A flow direction f ∈ F i determines the connection between the interior and exterior of each station by defining two sets of fence nodes: Those fence nodes V in f ⊆ V fn i allowing for flow into the station and those fence nodes V out f ⊆ V fn i allowing for flow out of the station.All other fence nodes do not allow any exchange of flow between the station's interior and exterior.
As an extension to the model given in [Hoppmann-Baum et al., 2021], we also decouple the pressure at those fence nodes without flow exchange.To do that, we adopt an idea initially mentioned in [Hoppmann-Baum, 2022]: We create a copy v of each fence node v, create a valve a in between the two nodes, and reconnect each artificial arc that was connected to the fence node v to its copy v instead.The construction is also visualized in Figure 4.Note that we add the copy v to V and the new valve a to A va .Thereby, we create all the corresponding variables and constraints.We further denote by a vaFn v the fence node We create a copy of the fence node and add a new valve between the fence node and its copy, the fence node valve of the node.The artificial station arcs that have been connected to the fence node are instead now connected to the copied node.
valve of fence node v, by A vaFn ⊆ A va the set of all fence node valves contained in the network, and by A vaOr = A va \ A vaFn the set of original non-fence-node valves of the network.
In addition to the flow restrictions, some fence nodes v have an upper pressure bound pexit v in case flow leaves the station over v for some flow direction, i.e. the fence node serves as an exit of the station.For all other fence nodes v, we set pexit v = pv .Finally, not all flow directions fit to all simple states.Hence, we are given the set Q ⊂ F × S of suitable combinations of flow directions and simple states.
For the model regarding flow directions and fence node valves, we demand for all network stations i ∈ I and all times t ∈ T that f ∈Fi x fd f,t = 1 (16) for all simple states s ∈ S i of all stations i ∈ I and all times t ∈ T that and finally, for all fence nodes v ∈ V fn i of all stations i ∈ I and all times t ∈ T that

Demand scenario
At each boundary node, we are given future requirements in terms of inflow and pressure, which we call the demands.A distinct set of demand values for all boundary nodes of the network is called a demand scenario of the problem.While for each future time step t ∈ T 0 , there is a required inflow value dv,t for each boundary node v ∈ V b , the prescribed pressure pv,t only exists for entries v ∈ V + with non-zero inflow values at time t.Furthermore, we do not have to meet them exactly but within a tolerance of ε = 1 bar for each value.However, there is no guarantee that meeting the demands is actually possible.For these cases, we allow a deviation from the demands, which we measure and penalize in the objective function.We call this deviation slack and capture it in the variables σ d+ v,t and σ d− v,t for the positive and negative inflow deviation and σ p+ v,t and σ p− v,t for the positive and negative pressure deviation for a boundary node v ∈ V b and time t ∈ T .When denoting by σd the maximum amount of inflow slack per boundary node and time step as well as by σp the corresponding maximum amount of pressure slack, we can formulate the model of the inflow demands for all boundary nodes v ∈ V b and a future time t ∈ T as and the model of the pressure demands for all entry nodes v ∈ V + with non-zero inflow demands dv,t and a future time t ∈ T as

Initial state of artificial elements
As stated at the end of Section 2.1, we are given an initial state for the network, specifying the pressures, flows, and modes of all nodes and elements of the original network.However, since the network station model is based on an artificial topology, there are no initial state values for the non-fence nodes inside the station, the artificial arcs, the station's simple state, and its flow direction.
For the inner station nodes, we do not need the initial pressure values as variables in the model.However, the initial pressures and temperatures are used as parameters in the artificial compressor model.For these values, we use the corresponding average values per quantity over the set of fence nodes.As the only exception, we use the initial quantities given for the original fence node for those nodes being fence node copies.Regarding the flow values, we do not use the initial flow values of artificial arcs.However, for the objective function, we need the initial flow on the fence node valve (see Section 2.10), which can be determined by using the flow balance equation (1) of the corresponding fence node.
Regarding the discrete states of the network stations as well as their contained artificial arcs, we aim to choose settings that fit the initial state of the fence nodes, fence node valves, and the compressor stations in the original network station topology.Since there might be multiple suitable settings, we only sort out non-fitting ones and include the final decision into the optimization model.
First, we determine a list of possible initial simple states S 0 .This is done by checking for each simple state of each network station if there would be a corresponding feasible initial state for the flow direction choice, as well as all the elements in the network station which fits to the initial pressure and flow values of the fence nodes and corresponding fence node valve.As a second step, we determine the set of initially running compressor units U 0 based on the active configurations of the compressors stations in the original network station topology.Finally, we determine the set C ar0 of those configurations of artificial compressors c, for which there are enough initially running compressor units contained in their set of usable compressor units U c , such that the configuration might be active.
Having these sets, we can define the model for the initial time step.From the above-defined constraints, we use (10a), ( 12), ( 14), and (15) also for t = 0.In addition, we add the following new constraints: x st s,0 = 0 ∀s ∈ S \ S 0 (21) a∈A arCp c∈C ar a :u∈Uc While the first two constraints ensure that only initially valid simple states and artificial configurations are chosen for the initial state, the last constraint demands for each initially running compressor unit u that there is at least one initially active artificial configuration c with u ∈ U c .

Objective function
The general goal of the model is to reduce the amount of changes necessary for fulfilling the given demands while trying to reduce unnecessary compressor unit usage.In contrast to the base model of [Hoppmann-Baum et al., 2021], where a multi-level model was formulated to ensure each solution has minimal slack values, we instead used a single-level objective function representing the sum of differently weighted terms.Another difference from the base model is that we do not perform a solution smoothing in a post-processing step but instead include this already in the model, which increases its complexity.
To build our objective function, we need to introduce additional variables to track the corresponding changes of the single entities and connect these with the quantities describing the network state.We start with the changes of valve modes.For each non-fence-node valve a ∈ A vaOr and time t 1 ∈ T , we capture a change in the binary variable δ va a,t1 and specify the corresponding model as δ va a,t1 ∈ {0, 1}, where t 0 = t 1 − 1 denotes the time step preceding t 1 .
For the regulator changes, we use the model introduced in [Hennings et al., 2021a] tracking changes of the regulator's mode and, in case it is active, the regulator's operation point in terms of changes in the incoming pressure, outgoing pressure, and flow.For a regulator a = ( , r) ∈ A rg and times t 1 ∈ T and t 0 = t 1 − 1, the model reads as ∀y ∈ {p , p r , q a } (25c) δ rg a,t1 ∈ {0, 1}, where δ rg a,t1 denotes the binary variable indicating a mode change at the regulator between times t 0 and t 1 , and δ rg-pl a,t1 , δ rg-pr a,t1 , and δ rg-qa a,t1 track the changes in the operation point for an active regulator for the incoming pressure, outgoing pressure, and flow through the element.
Next are the changes regarding the network stations and their corresponding elements.First, we introduce binary variables δ arc a,t , representing a change in the activity of the artificial arc a between time t and the previous time step.The corresponding constraints have the same structure as the valve constraints above and read for an artificial arc a ∈ A ar and times t 1 ∈ T and t 0 = t 1 − 1 as x δ arc a,t1 ∈ {0, 1}.For artificial compressors, we also count configuration changes as changes in the arc's activity.We only need to track if a configuration is newly activated, as turning it off either turns another configuration on or deactivates the whole compressor arc.Hence, we formulate for each configuration c ∈ C ar a of each artificial compressor a ∈ A arCp and times t 1 ∈ T and t 0 = t 1 − 1 the corresponding constraint as As in the base model, we track changes in the network stations' simple states but not in the flow directions.Since each station has exactly one simple state at each point in time, we capture the activation of a simple state s at time t in the binary variable δ st s,t and give the model for a simple state s ∈ S i of a station i ∈ I and times t 1 ∈ T and t 0 = t 1 − 1 as δ st s,t1 ∈ {0, 1}.
Finally, we capture the pressure changes at each fence node v between time t and the previous time point in the variable δ fn-p v,t and the flow changes at the corresponding fence node valve by the variable δ fn-q v,t .The corresponding model for a fence node v ∈ V fn i of network station i ∈ I and times t 1 ∈ T and t 0 = t 1 − 1 is given as Having defined all the change tracking variables, we are now able to formulate the objective function as a weighted sum of different penalty terms.Apart from the changes and slack values defined above, we also penalize the usage of compressor units by penalizing the usage of artificial configurations c weighted by their number of used compressor units u c .For each of the terms in the sum, we have specific weights defining the penalty costs: The demand slack cost per hour w σ-d , the pressure slack cost per hour w σ-p , the cost for each compressor unit running in an active configuration per hour w ru , the non-fence-node valve change cost w va , the regulator mode change cost w rg , the operation point change costs of active regulators for the incoming pressure w rg-pl , the outgoing pressure w rg-pr and the flow w rg-q , the artificial arc change costs w arc , the costs for activating a simple state w st s , which are specific values manually defined by the industry experts for each simple state s ∈ S, and the costs for changes of the fence node pressure w fn-p and flow changes on the corresponding fence node valves w fn-q .The complete objective is then given as .

Complete model
By combining all the above-defined variables and constraints with the objective function, we can formulate our gas flow problem P as the following MINLP model:

MIP-based heuristic algorithms
We aim at solving the transient gas flow problem presented in the previous section for challenging demand scenarios on large real-world-based instances.Since this is currently out of reach for black-box MINLP solvers, we create a specialized solution approach given as Algorithm 1.It is based on the idea of dealing with the two sources of complexity, the non-continuity of binary variables and the non-convexity of the constraints, separately: In the first step, we determine a promising and feasible set of solution values for all the binary variables, which we call a binary assignment of the problem.Afterwards, we transform the problem into an NLP by fixing these values and then solve this remaining problem with a black-box NLP solver.Finally, we return the overall best solution out of all successful NLP runs.If non of the runs found a feasible solution, the heuristic failed.Note that the algorithm could be easily parallelized by solving the NLP model immediately after finding a new binary assignment and in parallel to the ongoing search for binary assignments.Hence, solutions would be obtainable early during the execution, which is especially valuable in environments with strict time limit restrictions.Since we consider especially challenging and large instances, we purely focus on finding good feasible solutions fast, i.e., we do not aim for valid global lower bounds for the problem and thereby prove optimality of our solutions.However, note that in our approach, the black-box NLP solver usually provides, at least locally, optimal solutions for the corresponding sub-problem for a given binary assignment.
In the remainder of this section, we focus on finding binary assignments that yield high-quality solutions for the overall problem.Note that we sometimes claim that these binary-assignment-producing Algorithm 1: Base algorithm Data: An instance Φ of problem P Result: A valid solution for Φ 1 baSet ← findBinaryAssignments(Φ) 2 solSet ← solveNLP(baSet, Φ) 3 return chooseBestSol(solSet) algorithms solve the overall problem, which we use as a short form of explaining that we solve the problem by using them as function findBinaryAssignments in Algorithm 1.

Sequential mixed-integer programming
The core of our algorithm for finding binary assignments consists of a sequential algorithm solving linearized versions of the original MINLP problem, which we refer to as SMIP and which is given as Algorithm 2. After initializing a linearization for the non-convex Momentum equation (2b) in Line 2, we solve the remaining MIP problem.If successful, this yields a valid binary assignment in the sense that it respects all the constraints involving no continuous variables.We then check if the binary assignment was already found before and if the linearization error is small enough to fulfill the convergence criteria.If both are not the case and we have not reached the maximum number of iterations n it , we update our linearization based on the current solution and repeat the process by solving the updated MIP model.The algorithm finally returns the set of found binary assignments to be completed to full solutions for the original problem P. Note that given a valid binary assignment, the resulting NLP model might be infeasible and hence does not yield a feasible overall solution.

Linearization and convergence
As already mentioned above, the only non-linear term in the model is the friction term in the Momentum equation (2b) for pipes given as x,a,t + α p x,a,t • p x,t + α q x,a,t • q x,a,t .
To check if a solution of a linearized MIP model has converged against a non-linear solution, we use two tolerance parameters: The absolute convergence tolerance n tolAbs and the relative convergence tolerance n tolRel .Given these, we define that a solution with pressure values p and mass flow values q converged, if for all pipes and future time points the differences between the linearized friction function and the non-linear friction function with respect to p and q are smaller than the tolerances.Hence, we check in function isLinearizationErrorSmall used in Line 10 of Algorithm 2 if the following constraints are fulfilled: ,a,t + α p ,a,t • p ,t + α q ,a,t • q ,a,t + α 0 r,a,t + α p r,a,t • p r,t + α q r,a,t • q r,a,t Here, the value of the linearized friction function with respect to the variable values of the given solution is denoted by Λ lin and the corresponding non-linear friction function value by Λ orig .Furthermore, we can only check for the relative tolerance if the non-linear friction value is not zero.We ensure this by comparing the term against a corresponding epsilon value ε zero , which we define to be equal to 0.001 bar.
For the initial linearization determined in function applyInitialLinearization of Algorithm 2 in Line 2, we use the constant velocity approximation introduced by [Hennings, 2018], which was also used as initial linearization in [Hoppmann-Baum et al., 2021] and is given as , and q 0 := max{|q x,a,0 |, q min }.
Here, |v x,a,0 | denotes the absolute value of the gas velocity based on the initial state quantities using a minimal absolute flow value of q min to avoid α q x,a,t being zero.The advantage of this linearization is that it is symmetric in q and also zero for q = 0.
For all subsequent iterations of Algorithm 2, we determine the linearization based on the pressure values p and flow values q of a previous MIP solution by function linearizeBasedOnPreviousSolution in Line 12.As linearization, we use the tangential plane on the non-linear function in the point (p , q ), which is the first-order Taylor expansion at this point and is given as x,a,t • q x,a,t .

Comments on infeasibility
Algorithm 2 struggles with infeasibility in two different ways: First, the linearized MIP problems of each iteration might be infeasible even if the overall MINLP problem P is not.Second, the algorithm has no functionality to detect the overall infeasibility of the general problem P.
Despite those structural weaknesses, we still expect the algorithm to perform well on real-world-based instances, like those used in our computational experiments in Section 4. On these, infeasibility usually does not occur, as real-world gas transport networks are, in general, very flexible and allow for a broad range of different pressure, flow, and inflow values.In addition, the slack values can adjust problematic inflow and pressure demands.
The computational results we present in Section 4 confirm our expectations, as we found a feasible solution for each instance of the problem P and rarely encountered infeasibility when solving the linearized MIP models, see Table 5.

Reduced time horizon heuristics
To further accelerate the solving process, we introduce additional heuristics that find binary assignments by solving a series of models on time horizons of reduced size, where the size of a time horizon is defined as the number of future time steps.The heuristics are based on two well-known ideas for time-expanded problems: The rolling horizon and the aggregated horizon.While rolling horizon heuristics rely on iteratively solving the model for a subset of time steps that are gradually shifted to the future, the aggregated horizon combines certain time steps and afterwards completes the aggregated solution to one for the whole time horizon.
Note that we maintain the general sequential algorithm structure of Algorithm 2 for all of the heuristics below.

Rolling horizon heuristic
In our rolling horizon heuristic RH given in Algorithm 3, we only change a single line in contrast to the SMIP Algorithm 2: Instead of solving the linearized MIP for the whole time horizon in Line 4, we find a solution by using a rolling horizon approach of fixed length n rolHor .We start by building a reduced model in Line 8 that only considers the first n rolHor future time steps by not adding any variables or constraints for the remaining ones.After solving the reduced model, we save all variable values regarding the first of the considered future time steps for our overall solution in Line 12.We call this set of values the solution state of this time step and use it as the initially fixed state when building the reduced model in the next iteration.Note that we still use the overall initial time step 0 to determine the model's constant parameters, like, for example, the fixed temperature T a of a pipe a ∈ A pi .Finally, for the next iteration, we also include the earliest not already included future time step into our reduced time horizon, which thereby again has a length of n rolHor .We repeat this procedure until we reach the iteration containing the very last future time step t max .After saving all the solution states of this iteration in Line 14, we can create a solution for the entire time horizon by combining all the collected solution states in Line 18.
As a minor detail, we note that for all variables without a fixed initial state for time t = 0, for example, for the simple state variables x st , we use the model as defined in Section 2 for the first iteration and then add the there determined values for t = 0 to the corresponding solution state.For all subsequent iterations, we fix the values determined in the previous iteration for the initial state and can therefore ignore all the constraints added to the model for t = 0.

Aggregated horizon heuristic
For the aggregated horizon heuristic AH stated as Algorithm 4, we create a time horizon of reduced size n aggHor by combining certain time steps, resulting in longer but fewer time steps in the overall horizon.We represent this by a function ϑ : {0, 1, . . ., n aggHor } → T 0 = {0, 1, . . ., t max }, which maps the time steps of the aggregated time horizon to time steps of the original time horizon and is created as variable aggToOrigMap by calling the function mapAggregatedToOriginalTimeSteps in Line 1.A time interval between two subsequent time steps t and t + 1 in the aggregated horizon represents the combination of all original time intervals between the time points ϑ(t) and ϑ(t + 1).To combine a similar amount of original time steps for each time interval, we determine ϑ as c := t max n aggHor the larger number of original time intervals to combine c := t max n aggHor the smaller number of original time intervals to combine n := t max mod n aggHor the amount of larger aggregated time intervals and give an example as 0 2 4 6 7 .
Note that ϑ(0) = 0 and ϑ(n aggHor ) = t max hold.If n aggHor is a divisor of t max , each aggregated time interval combines c = c original time intervals.
Using the function ϑ, we can now derive a time aggregated problem instance Φ agg from the original problem instance Φ by calling the function createTimeAggregatedProblemInstance in Line 2 of Algorithm 4. Since all the elements sets of Φ and Φ agg coincide except of the time horizon, we are only missing a description of the time-dependent instance parameters.In a first step, we determine the inflow demands.As they represent the average inflow rate over a time interval, we determine the time aggregated inflow demands as the average of original inflow demands, weighted by the length of the corresponding time intervals, i.e., for the aggregated time t ∈ {1, . . ., n aggHor } and a boundary node v ∈ V b we define dv,t := For all the remaining time-dependent parameters, like, for example, the function τ (t) representing the time difference between t and the initial time step in seconds, the pressure demands at the entries, or the general pressure and flow bounds, we use for each aggregated time point t the parameter value of the corresponding original time point ϑ(t) in Φ and thereby complete the aggregated instance Φ agg .Using the SMIP Algorithm 2, we determine a set of valid binary assignments for the aggregated problem instance.Note that we use the parameter n itAgg to specify the maximum number of iterations to use in SMIP to solve Φ agg .In the following, we explain the steps to complete these to full binary assignments for the original problem instance Φ.The general idea is to fix certain variables in the model based on the aggregated binary assignment and thereby create the problem instance Φ fixed , which we then again solve by a variant of the SMIP Algorithm 2 to find the binary values for the rest of the original time horizon.Our set of variables to fix for creating Φ fixed is given in Line 4 of Algorithm 4.
Since we do not fix all the variables, the aggregated binary assignments created by solving Φ agg might be identical with respect to the variable fixations.In order to avoid unnecessary computation, we remove duplicated binary assignments by calling the function removeDuplicateBAsWithRespectToVarsToFix in Line 5 of Algorithm 4. Now, we create for each remaining binary assignment a variant Φ fixed of the original problem instance Φ by calling function createProblemInstanceWithPartiallyFixedBinaries in Line 8. When denoting the binary assignment by B and the value of a binary variable x in B by ς(x, B), we perform the following fixations: ∀t agg ∈ {0, 1, . . ., n aggHor } and variables x ∈ varTypesToFix : x ϑ(t agg ) = ς(x t agg , B)  // Return all full horizon binary assignment found be completing aggregated ones 11 return baSetFullHorizon In words, we do two types of fixations: First, we fix all binary variables of the corresponding types whose time steps are part of the aggregated time horizon to their corresponding solution values defined by B. Second and for all time steps in between those time steps from the aggregated horizon, we check if a variable has the same solution value with respect to B at the two surrounding time steps from the aggregated horizon.If it does, we also fix them to exactly that solution value.Otherwise, we do not fix its value.This means that we allow a variable value change between two aggregated time points to happen at any point in time in the corresponding original time interval.
After creating the instance Φ fixed , we again solve it by the SMIP Algorithm 2 using the parameter n itFull to specify its maximum number of iterations.However, we make two adjustments to the algorithm.First, we replace the applyInitialLinearization function with a newly defined function initializeLinearizationBasedOnAggSolution.In there, we still use the constant velocity approximation but determine the coefficient α q x,a,t for each time t, pipe a = ( , r) and end node x ∈ { , r} based on the pressures and flows of the aggregated solution that created the binary assignment.If t = ϑ(t agg ) for some t agg in the aggregated time horizon, we use the pressures and flows of that aggregated time step.Otherwise, we linearly interpolate the values.As a second adjustment, we provide a partial starting solution to each attempt to solve a linearized MIP, from which the solver tries to build a complete solution for the problem.We represent this adjustment by replacing the original solve function with the function solveMIPUsingPartialStartSolution.As the partial start solution, we provide the values for all binary variables x fd f , which determine the activity of a flow direction f , and m op a , which represent the mode for all fence node valves a.The values for each entity and time step are created based on the given binary assignment using the same algorithm as in function createProblemInstanceWithPartiallyFixedBinaries.Since both variables do not have a meaningful value for the initial time step, we use for all original time steps, which exist before the first time step with a corresponding aggregated time step, the variable value of that aggregated time step.
As a result of calling the adjusted SMIP in Line 9, we obtain a set of binary assignments for the original problem instance Φ.The union of these for all aggregated binary assignments is returned as the result of Algorithm 4. Note that the binary assignments obtained by completing different aggregated binary assignments are also different due to the guaranteed variable fixations on those original time steps with a corresponding aggregated time step.

Aggregated horizon heuristic using rolling horizon
As the final heuristic based on reduced time horizons, we propose a combination of the previous two algorithms and is given as Algorithm 5. We refer to it as ARH.
Algorithm 5: Sequential aggregated and rolling horizon heuristic algorithm (ARH) Data: An instance Φ of problem P Parameters: Size of the aggregated horizon n aggHor Size of the rolling horizon n rolHor Maximum number of iteration for the aggregated horizon MIP n itAgg Maximum number of iteration for the full horizon MIP n itFull Absolute convergence tolerance n tolAbs Relative convergence tolerance n tolRel Result: A list of valid binary assignments for Φ 1 ahUsingRh ← adjustAhByUsingRollingHorizonToSolveAggregatedMIP(n rolHor ) 2 return ahUsingRh(Φ,n aggHor , n itAgg , n itFull , n tolAbs , n tolRel ) The heuristic is nearly equivalent to the aggregated horizon heuristic given as Algorithm 4. The only difference is that we use in Line 3 the rolling horizon Algorithm 3 instead of the sequential mixedinteger programming Algorithm 2 to solve the problem instance Φ agg on the aggregated time horizon.We represent this in Algorithm 5 by the function adjustAhByUsingRollingHorizonToSolveAggregatedMIP, which returns a copy of Algorithm 4 adjusted by the mentioned change.This adjusted algorithm is then used to solve the given problem instance Φ.Note that we are given the parameter n rolHor , which specifies the length of the rolling horizon.

A dynamic node limit
When testing the above-given heuristics, we observed that during the solving process of the linearized MIP models, we find good primal solutions fast and then spend much time proving their optimality.This is typical behavior for branch-and-bound-based MIP solvers according to [Berthold et al., 2018].However, for our heuristics, we are only interested in the optimal (or at least a good) primal solution for each MIP.The corresponding dual bound is not relevant to us.Therefore, there is the potential to save time if we can end the solving process at the point at which the solver has found a good solution and changed its focus to improving the dual bound.
To implement the described behavior, we define a maximum limit L on the number of processed branch-and-bound nodes, which depends on the so far found feasible solutions and can be applied to the MIP solving processes in the heuristics.The general idea is to increase the node limit as long as we find new best primal solutions that significantly improve the primal bound.If this is no longer the case, we hit the dynamic node limit and stop the run with the current best solution.
At the start of the solving process, we set the limit to L = ∞.Whenever we find a new best feasible solution with objective function value O at node N , we update the node limit dynamically as follows: Let O last be the objective function value of the last feasible solution that updated the node limit and let O last = ∞ if no feasible solution has been found so far.Then update can be stated as In words, we update the node limit L whenever there is a significant relative improvement of at least n improve ∈ (0, 1) of the primal bound compared to its value at the last node limit update.The new limit has a minimal value of n limitMin and is otherwise set to n limitRel times the current number of processed branch-and-bound nodes.The minimal node limit avoids a premature termination of the solver and should cover the early stage of the solving process, in which the relative increments in the number of nodes are still relatively high.

Computational Experiments
The performance of the presented algorithms is tested by conducting a series of computational experiments.The corresponding set of instances is based on a real-world network with corresponding initial states as well as inflow and pressure patterns at the boundary nodes, which are provided by our project partner OGE [OGE].From this data we choose 40 especially challenging instances, by finding those with the highest amount of demand changes.After giving a description of the corresponding network features and describing the instance generation process, we specify the setup of the experiments as well as the used parameters.Then, we present the results of the two conducted experiments: First, we compare the results of the general sequential mixed-integer programming (SMIP) Algorithm 2 used in Algorithm 1 against those of a black-box global MINLP solver.Here, we focus purely on the solution quality to investigate if the algorithm is able to yield solutions close to the proven optimum.In a second experiment, we evaluate the general performance of the reduced time horizon heuristics presented in Section 3.2 while using the general SMIP approach as a reference.

Instance generation
The instances for testing our algorithms are based on a section of the real-world network of our project partner OGE [OGE], which was also used in [Hoppmann-Baum et al., 2021].It has been modified by replacing the original network station topology with the artificial model that was manually created by industry experts.Furthermore, the network outside of network stations is aggregated using the methods presented in [Lenz, 2021].We choose the demand scenarios defining our instances from a historical time period of 12 consecutive months.During this time, the network topology, the manually created station model, and the parameters of the applied network aggregation changed multiple times.Hence, we are given a slightly different network topology for each instance.An overview of the network characteristics is given in Table 1, in which the minimal and maximal amounts for each element type are given.Most noteworthy, the number of network stations changed from 7 to 10, as three regulators arrangements in the network were replaced by small stations.We note that compared to the instances used in [Hoppmann-Baum et al., 2021], a less aggressive network aggregation was used, resulting in more elements outside of network stations.1: Network topology statistics.Since the network and manually created station model change over time, we give for each quantity the corresponding minimal and maximal amounts.

|V| |V
In Table 2, we list parameters regarding the network stations.If at least one of these values changed over time, we give two sets of values for a station representing the corresponding minimal and maximal amounts per quantity.The network stations A to G are the same as in [Hoppmann-Baum et al., 2021], while the newly introduced stations replacing the regulators are named H, I, and J.
Artificial compressor model size Based on the specific network station characteristics in the given time period, we compare the actual sizes of the artificial compressors arc models used in the base model of [Hoppmann-Baum et al., 2021] and our newly introduced model from Section 2.5.4.We count the number of continuous and binary variables as well as the number of linear inequality constraints needed for each discrete time step.For the variables, we included the activation and flow variable per arc as well as all additional variables needed specifically for artificial compressors.The constraints we count for our model are those defined in Equations ( 10), (11), and ( 12), while we count all the constraints (31)-( 42) from [Hoppmann-Baum et al., 2021] for the base model.Since we count inequality constraints, equations are counted twice.The result is given in Table 3.We see that our model is smaller regarding all three categories, as it uses a slightly smaller number of binary variables and considerably fewer continuous variables and inequality constraints.Challenging demand scenarios To properly examine if the presented algorithms are suitable for application in time-critical industry environments, we search for particularly challenging instances.In [Hoppmann-Baum et al., 2020], where the gas flow model of [Hoppmann-Baum et al., 2021] was used to conduct a case study on hydrogen transport, the authors suggest that the difficulty of the problem corresponds to the number of element control changes needed to fulfill the given demand scenario.As the number of changes is not known a priori, we instead look for demand scenarios featuring large demand changes for single boundary nodes since we expect these changes to induce inevitable changes in some network elements' control.From a cooperation with our project partner OGE, we attain historical network demands in hourly resolution over a time period of 12 months, covering the majority of the year 2020.The time period is nearly consecutive, having only three periods of 2, 6, and 11 missing days due to problems during the data creation.
To find demand scenarios with large demand changes from within this time period, we first specify the desired time horizon length for our instances to be 12 and 24 hours.For both, we then enumerate all time horizons of the corresponding length and sum the differences between the initial and end time point of the two quantities used as future demands in Section 2.8: The boundary node inflow values, which are given as norm volumetric flow Q 0 in 1000m 3 /h, and the pressure in bar for entries with positive inflow.The norm volumetric flow is defined as Q 0 = q/ρ 0 , where ρ 0 is the norm density depending on the gas mixture.We then create a final score ∆ for the amount of changes of a time horizon by adding the two sums, using the quantities' values for the given units, and weighting the pressure changes by a factor of 100.Finally, we choose from each of the two lists of time horizons a set of 20 demand scenarios by iteratively selecting from the remaining time horizons the one with the highest ∆ value, which does not overlap with any of the already selected ones.Overlapping is hereby defined as sharing at least one full hour in the time horizon.If one time horizon ends at time point t and a second one starts at t, they do not overlap.
An overview of the ∆ values for each of the 20 demands scenarios for each time horizon is given in Figure 5.The largest found score value is about 6600, which would be equivalent to a total amount of inflow changes of 6600 1000m 3 /h or a total amount of pressure changes of 66 bar at entry nodes over the whole network.The 40 demand scenarios, together with the corresponding network topology and initial network state at the start of the time horizon, define the 40 instances we use as our test set for the computational experiments.When using a time horizon whose time steps are larger than 1 hour, we aggregate the hourly demands in the same way as in the function createTimeAggregatedProblemInstance used in the aggregated horizon heuristic presented in Algorithm 4.

Setup and parameters
In this section, we give an overview of the computational hardware setup and solver choice, the applied model adjustments to improve its numerical properties, and specify the parameter used for the computational experiments.

Computational setup
The computations were executed on a cluster using 4 cores and 40 GB of RAM of a machine composed of two Intel Xeon Gold 5122 running at 3.60 GHz.The MIP problems are solved by Gurobi in version 9.5.1 [Gurobi Optimization, LLC., 2021], for which we change the general solution strategy to focus on finding primal solutions by setting the MIPFocus parameter to 1.As the NLP solver, we used Ipopt in version 3.14.3[Wächter and Biegler, 2006].Finally, we used Scip in version 8.0.0 [Bestuzheva et al., 2021] as the global MINLP solver.For Scip, we change the feasibility tolerance for constraints numerics/feastol to 10 −3 , which is even higher than the corresponding default value of 10 −4 for Ipopt, to improve the acceptance rate for start solutions.We accessed both Ipopt and Scip via Gams in version 38.1.0[GAMS Development Corporation, 2022].

Coefficient scaling and rounding
To improve the numeric properties of the models, we use similar ranges for the continuous variables by scaling them such that pressure variables are used in bar and flow values in kg/s.Furthermore, we apply the following modification to the pipes' Continuity equation (2a), their Momentum equation (2b) in the original or the linearized variant, and the linearized maximum power constraint of artificial configurations (11c): First, we set all coefficients to zero, whose absolute value divided by the largest absolute value coefficient is smaller than 10 7 .Afterwards, we scale all coefficients of the corresponding constraint such that the smallest non-zero coefficient has an absolute value of 1.

Parameters
In the following, we list the parameters used in the model, in the heuristic algorithms, and as resources limits for the different solvers.
Model parameters For the model, we specify the maximum slack parameters to be 2 bar for the pressure slack.In the case of the inflow, we restrict the possible maximal deviation from the demands to 50%.Hence, we use, instead of a fixed upper bound σd on the inflow slack, an upper bound specific to each boundary node and time step of σd v,t = 0.5 • | dv,t | for all boundary nodes v ∈ V b and future time steps t ∈ T .Furthermore, we set the minimum flow value used to determine the absolute velocity in the initial linearization of the friction term to q min = 100 • 1000m 3 /h.
Note that the multipliers for the flow change and the hourly flow slack are given with respect to norm volumetric flow Q 0 in 1000m 3 /h.Hence, we multiply them by a factor of 3.6/ρ 0 to get the corresponding mass flow equivalent.The costs for turning on a simple state of a network station are manually created by the industry experts at OGE, which is why we can only give the range of values for them.
Algorithm parameters Each heuristic algorithm defined in Section 3 is given a set of parameters to adjust them.First, we specify the converge tolerances for the friction term linearization as n tolRel = 0.001 n tolAbs = 0.01 bar.
Next, the maximum number of sequential mixed-integer programming iterations we use is Finally, we do not fix the parameters defining the reduced time horizon size globally but define different variants of the heuristics algorithms based on these.Hence, we denote the rolling horizon heuristic of Algorithm 3 with a rolling horizon size of n rolHor = X by RXH, the aggregated horizon heuristic of Algorithm 4 with aggregated horizon size of n aggHor = X by AXH, and the aggregated and rolling horizon heuristic of Algorithm 5 using an aggregated horizon of length n aggHor = X and a rolling horizon of length n rolHor = Y by AXRYH.
Resource limits As final parameters, we define the resource limits of the different solving processes, which are given in Table 4.There we distinguish five different model variants: MIP with the dynamic node limit defined in Section 3.3, MIP without this node limit, an NLP run to complete binary assignments, a fast NLP variant having a time limit of 5 minutes, and finally, the MINLP run.The given parameters for the dynamic node limit specify that the smallest node limit to stop is 1000 and that a significant primal solution upgrade featuring an objective function value improvement of at least 2% increases the node limit to 150% of the current number of processed nodes.Note that all variants are additionally restricted by the general memory limit of 40 GB.

Solution quality of general SMIP
In our first experiment, we check the quality of solutions obtained by using the general sequential mixedinteger programming SMIP Algorithm 2 to determine binary assignments in Algorithm 1 and complete these to a full MINLP solution by an NLP solver.The lower bound to compare the solutions against is obtained via a global MINLP solver.Due to the size of our network and the general complexity of the Variant Rel.Gap Time limit n limitMin n limitRel n improve MIP w/node limit 0.0 ∞ 1000 1.5 0.02 MIP wo/node limit 10 Table 4: Optimality conditions and resource limits for the different model variants.We distinguish solving a MIP with and without the dynamic node limit defined in Section 3.3, solving an NLP with and without a time limit, and solving an MINLP.Since the overall run time of the solving process has no priority for this experiment, we use for the single MIP models the parameter variant without the heuristic node limit, and for the NLP the general variant without a time limit, see Table 4. Also, we use the solution obtained from the SMIP algorithm as a starting solution for the MINLP.The result is given in Figure 6.First, we notice that for only 12 out of the 40 instances, the MINLP was solved to optimality.In 7 cases, the time limit of roughly one week was hit, while for the remaining 21 instances, the memory limit of 40 GB was exceeded.As a consequence, we do not always have the best possible lower bound available.However, the average gap of 3.2% of the SMIP solutions to that lower bound is still rather small and even closer to the average gap of 2.3% of the MINLP itself.Using the best found primal solutions of the MINLP as a reference for the highest possible lower bound values, the average gap for the SMIP solutions is 0.9%.Regarding run time, Algorithm 1 using the SMIP always finishes in less than 30 seconds, while the geometric mean run time of the MINLP is larger than a day.In summary, the SMIP solutions are very close to the global optimum and can be computed in only a fraction of the global solver's execution time.

Performance of reduced time horizon heuristics
In our main experiment of the paper, we evaluate our proposed algorithm of using a reduced-time-horizon heuristic of Section 3.2 for creating binary assignments to be completed in Algorithm 1.As a reference approach, we use the SMIP for the binary assignment creation, which was shown to find close-to-optimal solutions in the previous section.Not having to compare against a global MINLP solver enables us to solve instances with a time horizon consisting of 12 equally large time steps of 1-hour length for the 12hour horizon and 2-hour length for the 24-hour horizon.Depending on the actual instance, the MINLP model, which is never formulated and tackled as a whole model here, has the following sizes (rounded to 2 significant digits): After conducting some preliminary experiments, we determined for each of the reduced horizon heuristics the smallest time horizon lengths, which still give reasonably good results, and the largest time horizon lengths, at which the run times are still reasonably short.As a result, we choose the following variants to test for each of the heuristics: The rolling horizon heuristics from R1H to R4H, the aggregated time horizon heuristics A2H, A3H, A4H, and A6H, and for the last two aggregated horizon lengths also the heuristic using the rolling horizon heuristic for some of the sub-models with horizon lengths of 1 to 3, resulting in the variants A4R1H, A4R2H, A4R3H, A6R1H, A6R2H, and A6R3H.For all variants, we solve the MIP models using the dynamic branch-and-bound node limit setting and solve the NLP using a time limit, see Table 4.For the reference SMIP approach, we ensure the highest possible solution quality by not using the node limit for MIP or the time limit for NLP.
# Table 5: Statistics for the used heuristic variants.The first two columns state the name of the heuristic and the number of instances for which it fails to find a feasible solution.The next four columns state the average number of MIP models solved per instance, followed by the corresponding average number of those MIPs finished with an optimal solution, finished due to reaching the dynamic node limit, or finished with proven infeasibility of the corresponding model.Next, we have a similar block of four columns for the average number of NLP models solved per instance.First, the total average is stated (which is equal to the average number of binary assignments found), then the average number of converged NLPs, the average number of NLP finishing with a feasible solution at the given time limit, and finally, the average number of failed attempts to find a feasible solution within the time limit.The final two columns depict the geometric mean run time for the NLP in seconds, first for a single NLP run and second for the summed run time of all NLP runs per instance.If a heuristic did not produce any binary assignments for an instance, the NLP run time was set to 1 second.
We present the results averaged over all instances for each heuristic variant in two pairs of figures.Each figure shows the position of each variant with respect to the run time and the average gap to the virtual best solution.In the Figures 7, the run time is determined as the geometric mean run time of all   the instances and in the Figures 8, we use the maximum run time.If a heuristic failed to find a solution for a specific instance, we used a gap value of 100%.In Table 5, we gave an overview of the number of fails as well as the average results for the single MIP solves and NLP solves per heuristic variant.In the figures, we connected those heuristics being Pareto efficient by a black line, where Pareto efficient heuristics are those not having another heuristic with a smaller average gap and, at the same time, a shorter run time value.
In addition to the individual heuristics displayed in the upper Figures 7a and 8a, we also present in the lower Figures 7b and 8b results for the combination of two different heuristic variants.For each combination and instance, we use the best solution found in any of the two heuristics and add up the overall run time of both approaches.As a reference, we also included all Pareto efficient individual heuristics from the corresponding upper pictures.For the combination, we always choose a pair of a pure rolling horizon heuristic and a heuristic using an aggregated time horizon.We do this since we expect the two fundamentally different approaches to have a higher probability complementing each each other in finding high-quality solutions.For the combination, we chose all four rolling horizon heuristics.For the time-aggregation-based heuristics, we used those that are part of the Pareto front for the geometric mean and the maximum run time figures, i.e., A6H, A6R1H, A6R2H, and A6R3H.Note that the combined heuristics in the figures have the same color for using the same rolling horizon heuristics and the same symbol for using the same time-aggregation-based heuristic.
The figures show that the reference SMIP approach is the overall best approach in terms of solution quality with an average gap of 0.0%, as it yields the best solution for all of the instances.However, it is much slower than every of the reduce time horizon heuristics both on average and regarding the maximum run time.This is mainly due to the fact that, on average, 1.5 out of the at most 5 solved MIP models reached the 1-day time limit before proofing optimality.For the reduced time horizon heuristics, the solution quality is, in general, very good: The average gap to the virtual best solution is below 20% for all considered variants and even below 10% for all but two.The best average gap apart from the general SMIP approach was achieved by the A6H heuristic with a gap of 1.07% and 17/40 instances, for which the best solution has been found.Two very close variants are A6R2H and A6R3H, which are close in terms of the average gap but have better run time values.Hence, the summary for the individual variants is to choose an aggregated horizon heuristic with a long time horizon and either solve it directly for the best results or internally use the rolling horizon approach for shorter run times at the cost of a slightly worse solution quality.When comparing the geometric mean and the maximum run time, the pictures are very similar, with only two heuristic variants being in only one of the corresponding Pareto fronts.The maximum run time values are between 2 to 5 times higher for all the heuristics except A2H, where the factor is over 7.For the two variants A6R2H and A6R3H with the lowest run times of those with an average gap below 3%, the maximum run time is about half an hour, which is very fast for these types of problems.
Using two heuristics in combination improves the overall solution quality as expected, with 6 of the 16 combinations having an average gap better than the best individual solution and all of the combinations involving A6H reaching an average gap below 1%.However, this comes at the cost of rather large run time increases.Only the fastest of the combined heuristics has a better geometric run time than the one of the slowest Pareto efficient reduced time horizon heuristic.However, in that larger run time segment, several Pareto efficient combinations of the A6R3H and A6H heuristic combined with the rolling horizon heuristics can be found for both run time metrics.The best of them decrease the gap to the SMIP variant to a minimum of 0.42% for the two combinations of R3H and R4H with A6H while still being much faster than SMIP itself.
The reason for the overall bad run time for the combination of two heuristics is the rather long solve time for the NLP model, as displayed in the last two columns of Table 5.For heuristics like R1H and R2H that determine binary assignments very fast, the geometric mean NLP time per instance is equal to 96.5% and 71.5% of the geometric mean overall run time, respectively.Hence, testing more binary assignments for yielding good overall solutions to the problem is apparently a considerable time investment.
To compare our algorithm against a procedure similar to the one proposed in [Hoppmann-Baum et al., 2021], we also tested the SMIP approach with only a single iteration, i.e., n it = 1.This is close to the algorithm proposed there, as the linearization used for the model determining the binary solution values is only updated in the rare case of failure in the subsequent procedure searching for a non-linear solution for the found binary assignment.The SMIP heuristic with n it = 1 has an average gap of 42.89% to the virtual best solution and is, therefore, worse than all of our proposed heuristic variants.For 14 out of the overall 40 instances, it reached the time limit of 1 day.This leads to a geometric mean run time of 9541.0 seconds, which is again considerably larger than those of all our heuristics variants.Hence, we can conclude that our proposed algorithms are clearly superior.

Outlook
In this article, we presented several improvements for solving the control optimization problem on transient gas networks in a time-critical environment.First, we introduced a new model for the artificial arcs representing the compressors, which is configuration-based and closer to the original compression capabilities while using fewer variables and constraints.Furthermore, we propose a new algorithm that finds better solutions faster.It first determines a set of solution values for all binary variables, which we call a binary assignment, and afterwards solves a corresponding non-linear program to complete it to a full solution of the problem.The binary assignments are obtained based on sequential mixed-integer programming, two heuristics based on reduced time horizons, and a specialized dynamic node limit.Our approach was tested on particularly challenging real-world instances featuring a large amount of supply and demand changes.The results confirm that we are able to find solutions close to the global optimum while being very fast compared to black-box MINLP solvers and comparable approaches from the literature.
We see two main reasons for the success of our heuristic approach for the given problem.First, we recognized that from the multitude of possible future control recommendations, only a few combinations fit a given scenario.This is especially true for the given objective function that mainly punishes control modes and network state changes.Second, the network's pipes act as gas storage and can therefore compensate for a local shortage or surplus of gas for a short amount of time.For this reason, changing the control of an element a little too early or too late is often feasible and only causes a small objective function value penalty.This helps for both of the heuristics: The rolling horizon heuristic results in too late control changes due to not having a full view of the future.In the case of the aggregated horizon, the time points in which control decisions can be made are limited.Hence, they may be a little late or a little early.Due to the pipe storage, solutions with these control time offsets are likely to be similar to one of the few reasonable control recommendations while only having a slightly worse objective function value.
Future research on the overall topic can continue in many directions.Model-wise, we could further refine the compressor model to guarantee that the corresponding feasible region is an actual relaxation of the original feasible region.Another possibility to improve the accuracy would be to also update the linearization of the power bound during the sequential mixed-integer programming and include the original non-linear formula in the final non-linear programming model.However, as we have seen in the results, the run time for solving the NLP models is already quite high.Therefore, it might be beneficial to just complete the most promising binary assignments to full solutions, for example, by comparing their corresponding linearized MIP solution value, or to only complete one of a set of similar binary assignments.Regarding the instance set used in our computational experiments, it would be interesting to examine if heterogeneous step lengths in the time horizon lead to more complex solutions or to evaluate our algorithm in the context of pure hydrogen networks similar to those presented in [Hoppmann-Baum et al., 2020].

Variable
Meaning Unit p v,t ∈ R ≥0 pressure at node v ∈ V bar q v,a,t ∈ R flow into or out of pipe a = (l, r) ∈ A pi at end node v ∈ {l, r} kg/s q a,t ∈ R flow over arc a ∈ A \ A pi kg/s Figure2: Relative error in compression power.We compare the power P sum needed for a serial compression using two separate compression processes with ratios ψ 1 and ψ 2 to the power P one needed by our approximation of using a single compression process with ratio ψ 1 • ψ 2 .We plotted the relative error (P one − P sum )/P sum for compression ratios ψ 1 and ψ 2 up to 2.0. Figure created with [Plotly Technologies Inc., 2015].

Figure 3 :Figure 4 :
Figure 3: Feasible region comparison of the artificial compressor base model against our new one at the example of one artificial compressor arc.As a reference solution, we plotted the feasible region polytopes of the configurations of the associated compressor station from the original network topology in red.The polytopes for the artificial compressor arc are plotted in blue, where Figures (a) and (c) show the base model polytopes and Figures (b) and (d) show the new model polytopes.For the new model, each polytope represents one artificial configuration to choose from.Since the base model does not feature configurations, we plotted one polytope for each feasible compressor unit combination on the artificial arc.The pressure bound interval used for both end nodes is [20 bar, 80 bar].For the constant incoming pressure needed for the base model, we use the median value of 50 bar.For a better impression of the three-dimensional polytopes, we show them from two different camera positions in View 1, visible in Figures (a) and (b), and View 2, visible in Figures (c) and (d).Figures created with [Plotly Technologies Inc., 2015].

Figure 5 :
Figure 5: Change scores of the 20 non-overlapping demand scenarios with the highest change score ∆ for each of the two lengths of the time horizon.Each bar represents one selected instance and is divided according to the pressure and flow share of ∆, where the pressure is scaled by 100.The instances are sorted by change score.Figures created with [Plotly Technologies Inc., 2015].

Figure 6 :
Figure 6: Comparison of the gap of the solutions created by the SMIP Algorithm 2 used in Algorithm 1 with respect to the lower bound found in the MINLP run and the gap of the MINLP itself.Each pair of bars represents one instance, sorted by the gap of the SMIP solution.The values are given in percent.For the instances marked by a star, the MINLP proved optimality.For those instances marked by an hourglass, the MINLP time limit was hit, while for instances marked by a black triangle pointing upwards to a horizontal bar, the MINLP hit the 40 GB memory limit.Figure created with [Plotly Technologies Inc., 2015].
Combine heuristics and Pareto points of the individual heuristics

Figure 7 :
Figure 7: Position of the heuristic variants in a geometric-mean-run-time vs. gap-to-virtual-best-solution plot.If a heuristic does not find a solution for some instance, we count the corresponding gap as 100%.The Pareto efficient points are linked by a black line, forming the Pareto front.Figures created with [Plotly Technologies Inc., 2015].
Figure 8: Position of the heuristic variants in a maximum-run-time vs. gap-to-virtual-best-solution plot.If a heuristic does not find a solution for some instance, we count the corresponding gap as 100%.The Pareto efficient points are linked by a black line, forming the Pareto front.Figures created with [Plotly Technologies Inc., 2015].

Table 2 :
Network station statistics.Since the network and manually created station model may change over time, we give for each station and each quantity the corresponding minimal and maximal counts.If there are no changes in a station, we only give one set of values.

Table 3 :
Artificial compressor model sizes of the base and our new model based on the used instance set.Numbers vary over time.The corresponding minimal and maximal numbers originate from the same point in the overall time period.
,t ∈ R inflow into boundary node v ∈ V b kg/s m op a,t ∈ {0, 1} selection of open mode for a ∈ A va ∪ A rg 1 m cl a,t ∈ {0, 1} selection of closed mode for a ∈ A rg 1 m ac a,t ∈ {0, 1} selection of active mode for a ∈ A rg 1x arc a,t ∈ {0, 1} activity of artificial arc a ∈ A ar 1x fwd a,t ∈ {0, 1} selection of forward direction for a ∈ A arBiRg 1x bwd a,t ∈ {0, 1} selection of backward direction for a ∈ A arBiRg 1 ≥0 positive pressure slack for boundary node v ∈ V b barσ p− v,t ∈ R ≥0 negative pressure slack for boundary node v ∈ V b bar ≥0 change of pressure of fence node v ∈ V fn i for i ∈ I bar δ fn-q a,t ∈ R ≥0 change of flow over valve of fence node v ∈ V fn i for i ∈ I kg/s Table A.1:List of all used variables, specifying their domain, meaning and unit.Note that 1 bar = 10 5 Pa.