Model order reduction for gas and energy networks

To counter the volatile nature of renewable energy sources, gas networks take a vital role. But, to ensure fulfillment of contracts under these circumstances, a vast number of possible scenarios, incorporating uncertain supply and demand, has to be simulated ahead of time. This many-query gas network simulation task can be accelerated by model reduction, yet, large-scale, nonlinear, parametric, hyperbolic partial differential(-algebraic) equation systems, modeling natural gas transport, are a challenging application for model order reduction algorithms. For this industrial application, we bring together the scientific computing topics of: mathematical modeling of gas transport networks, numerical simulation of hyperbolic partial differential equation, and parametric model reduction for nonlinear systems. This research resulted in the morgen (Model Order Reduction for Gas and Energy Networks) software platform, which enables modular testing of various combinations of models, solvers, and model reduction methods. In this work we present the theoretical background on systemic modeling and structured, data-driven, system-theoretic model reduction for gas networks, as well as the implementation of morgen and associated numerical experiments testing model reduction adapted to gas network models.

is (and was [138]) never sufficient. To this end we evaluate customized model reduction techniques for an established class of gas network models.
This work and the associated software platform are an effort to determine which model reduction methods are suitable for enabling digital twins [58,87,121] of gas networks. Depending on the mathematical model and quantities of interest, the twin may contain redundant or superfluous information with respect to the simulations. Therefore, model reduction compresses the twin to a matched surrogate model, which is sufficiently accurate in the chosen operating region.
The swift numerical simulation of gas network twins by reduced order modeling is highly relevant, not only due to the transition towards renewables at the time of writing, which is underlined by the research projects MathEnergy 1 (Mathematical Key Technologies for Evolving Energy Grids) [29] that the authors are part of, and TRR154 2 (Mathematical Modelling, Simulation and Optimization using the Example of Gas Networks) [85], but also because of the intriguing numerical problem of model reduction for hyperbolic, nonlinear, coupled, parametric, multiscale partial differential-algebraic equation systems.
If relevant intraday demand changes occur, established steady / stationary / static simulations may not be sufficient anymore [42]. The basic model for the simulation of unsteady / dynamic / transient flow processes in gas pipelines is based on the onedimensional (isothermal) Euler equations, originally introduced in [55], and popularized in [104] as well as in [80,81] around the same time. A practical extension in the context of gas networks is the repetitive modeling approach [37], which enables a modular construction. For extensive details on gas network modeling, see the works [19,33,45,94,105,123,128,131], and for a concise summary of the overall approach we recommend [7]. Furthermore, a system-theoretic approach to gas networks is discussed in [35,54], and results on boundary reachability (controllability) and observability for this class of models have been derived in [8,9].
In this work, we conceptually combine these previous approaches, by using systemtheoretic but data-driven methods that are structure-preserving. The utilized data-driven assembly of the system-theoretic operators, central to the employed methods, is also a partial answer to the challenges posed in [57,Remark 5.10]; while structure preservation means in this context, retaining (particularly not mixing) the discretized physical variables in the reduced order model. Furthermore, we note that from this work's point of view, [130,144,145] are concerned rather with (valuable) model simplifications than model reduction.
To avoid the analytically most complex aspects of the gas network model -the nonlinearities -one could linearize the model around an operating point. Yet, the different nonlinearities (i.e. friction, compressibility and compressors) are unlikely to have a compatible operating point for a wide range of scenarios. Furthermore, linearized and simplified N q Dimension of mass-flux space s q Supply node mass-flux n p Dimension of reduced pressure space d p Demand node pressure n q Dimension of reduced mass-flux space models of gas flow have limitations with simulations of real scenarios [61], [134,Ch. 7]; hence, we use a nonlinear model. Since there is no general theory for model reduction of nonlinear systems, and a high degree of modularity in the gas network modeling process, model reduction algorithms have to be compared heuristically to determine their applicability. As a result of this reasoning and a demand for gas network simulation software tools [32,69], a platform named morgen (model order reduction for gas and energy networks) was designed with the goal to compare different models, solvers and reductors. The morgen 3 platform is an open-source project, and designed in a configurable, modular, and extensible manner, so that modeling, discretization or model reduction specialists can compose and compare their methods fairly. In summary, this work contributes a full, but also fully modular, modeling, model reduction and simulation open-source software stack for gas networks, and potentially other energy network systems (i.e. district heating networks, water networks), which brings together research results from various disciplines.
Overall, this work is organized as follows: In Sect. 2 the gas network model, simplifications, non-pipe elements, a relation to port-Hamiltonian theory, and obtaining a steadystate initial condition are described. Section 3 and Sect. 4 outline the general model reduction idea and propose five reduction method classes. The design and features of the morgen platform are summarized in Sect. 5, followed by three sets of numerical experiments in Sect. 6. We conclude by an outlook (Sect. 7) and an evaluation of our findings in Sect. 8. A list of recurring symbols is found in Table 1.

The transient gas network model
The goal of this section is to describe the partial differential-algebraic equation model of a gas network as an input-output system that maps boundary values to quantities of interest. First, the model for a single pipeline is summarized, which is based on the isothermal Euler equations of gas dynamics [104,106]. Then, it is generalized to a network of pipes, and simplified compressors are added. Additionally, a connection to energy-based modeling is made.
Even though further non-pipe elements are common in gas networks, such as resistors, coolers, heaters, valves and control valves [45,96], we prioritized compressors to focus on the model reduction aspect on a macro scale. Moreover, the practical numerical problems of scale homogenization, spatial discretization, index reduction and steady-state approximation are discussed in this section.

The gas pipeline model
The principal building blocks of gas transport networks are pipelines or ducts. Since the length of pipes exceeds their diameter by far (L > 500d, [81]), a spatial one-dimensional model suffices. We model gas flow in a (cylindrical) pipe of length L connecting two junctions by the isothermal Euler equations: which determine the evolution of the coupled pressure p(x, t) and mass-flux q(x, t) variables. The physical dimension of the pipe enters as its diameter d and the derived crosssection area S = π 4 d 2 , which is assumed constant, ignoring the influence of temperature and pressure on the pipe walls. These coupled partial differential equations (PDE) can also be characterized as a nonlinear, two-dimensional, first-order hyperbolic system of conservation laws: the pressure p preserves continuity, while the mass-flux q conserves momentum.
Following [64,79,108] and [19,Sect. 2.1], the inertia term has been neglected due to a low Mach number m 1, which leads to the ISO2 model in the TRR154 classification [33], also known as friction-dominated model [24,Sect. 3.2.1]. Furthermore, we assume a turbulent flow with a Reynolds number exceeding Re 10 5 [40,55], neither line breaks or valve closings happen intraday (to preclude associated shocks [34]), and low-frequency boundary values [10,108], which in this work are (but generally not limited to) the supply pressure and demand mass-flux, due to frequent use in literature, and use-cases like guaranteed demand pressures [22,63].
In (1), the linear reaction term describes the effect of gravity (with standard gravity g ≡ 9.80665 m s 2 ) due to the pipe height h(x), while the nonlinear reaction term models loss of momentum due to friction at the pipe walls, specified by the (Darcy-Weisbach) friction factor λ 0 := λ(d, k, Re 0 ), given a pipe roughness k, and an estimated mean Reynolds number Re 0 , see [19,Sect. 2.2]. 4 This friction term is principal to the accuracy of the gas pipeline model [33,61,93,107].
In this model variant, a (globally) constant mean compressibility factor z 0 := z(p 0 , T 0 ) ∈ R is assumed [35,62,106,111], as well as a constant gas state γ 0 := R S T 0 , whereas the temperature T 0 and the specific gas constant R S are treated as parameters (see Sect. 2.6). To this end, the steady-state pressurep =: p 0 is used to compute z 0 , via heuristic formulas based on the Virial expansion [25], [19,Sect. 2.3]. 5

Homogenizing scales
The SI-based units for pressure and mass-flux are [Pa] and [kg/s], respectively. This introduces a difference in scales of five orders of magnitude between the variables p and q, and hence induces numerical problems. To counter this multiscale structure, we simply rescale the pressure from [Pa] to [bar] which conveniently comprises a factor of 10 5 . Nonetheless, the model still consists of two interacting physical variables, hence the model still has to be treated as a coupled system, however, without numerical multiscale issues.

The gas network model
Given the model for a single pipe from the previous section, a (gas) network of pipes can be encoded by a finite directed graph, which is a tuple G = (N , E) of finite sets symbolizing nodes N , and oriented edges E. The edges correspond to pipes, while the nodes represent the junctions connecting pipes. The connectivity of the network is the relationship between edges and junctions, given by the incidence matrix A ∈ {-1, 0, 1} |N |×|E| , a map from edges to nodes, such that: Note, that the orientation of the edges is not enforcing the dynamic flow direction of the gas, but is necessary to determine the complexity and boundary of the overall networked model [19,53].
We introduce the notation |A| for the component-wise absolute value of a matrix. Using this absolute value, the following partial incidence matrices associating edges entering and leaving nodes respectively are defined similar to a Heaviside function: Next, based on this connectivity, certain conservation properties are enforced to maintain a network balance, and thus ensure physical relevance of the gas network model. Specifically, the Kirchhoff laws are applied to the network in vectorized (or rather matricized) form [18,132]: 1. The sum of in-and outflows (mass-flux) at every node (junction) is zero: This means that no gas gets lost in transport from one pipe to the next, with the exception of boundary nodes. Hence, a vector of flows q ∈ R |E| applied to the incidence matrix equals the (out-)flow at the boundary (discharge) nodes d q : with N D ⊂ N denoting the subset of boundary nodes, which only connect from one node respectively, but not to any node. 2. The sum of directed pressure drops in every fundamental loop is zero: Fortunately, an equivalent representation [143,Ch. 7.3] can be used, which resolves implicitly. It remains to ensure that the nodal pressures at the in-flow boundary (supply) nodes are associated to the boundary function s p : R → R |N S | , which are mapped to the network via B s ∈ {0, 1} |N S |×|E| : with N S ⊂ N denoting the subset of boundary nodes, which only connect to one node respectively, but not from any node, and the reduced incidence matrix with all rows associated to supply nodes removed. Given a connected and directed graph representing a gas network topology, with the dynamic flow in the pipe edges modeled by the PDE (1), then yields a partial differentialalgebraic equation (PDAE) due to the above constraints.

Discretization and index reduction
Next, we delineate the discretization of the spatial differential operators and reduction of the (P)DAE index in the networked system, yielding a system of ordinary differential equations (ODE). Eventually, the remaining discretization of the temporal differential operators is addressed.
We explicitly do not use the decoupling approaches from [12] or [11], as the former employs linearization and hence does not fit this setting, while compared to the latter, our equivalent analytic index reduction is more convenient here.
The partial differential(-algebraic) equation is discretized using the method of lines: First in space, then in time, yielding a (nonlinear) dynamic system. For the spatial discretization a first-order upwind finite difference scheme is utilized [10,133]. We select (only) two points for each of the k pipes with length L k , namely the start (· R ) and end point (· L ): The matter of short, long and varying lengths L k is addressed in Sect. 2.4.3.
For each pipe, this leads to the following equations: Now, different choices for (· * ) are surmisable. Subsequently, two specific combinations of p * and q * will be discussed: The midpoint discretization [18,50,51,53,145], and the leftright discretization [48,52,115] resulting in (implicit) ODEs. For an error analysis of these two discretization variants, see [128]. In the following, we describe a unified approach of deriving these index-reducible discretizations. For notational ease in the coming subsections, a vectorized form for the above (networked) system including its constraints is given by: using, thus resolving, the constraint A T 0 p + B T s s p = p Rp L , as well as d 0 := 1 γ 0 z 0 ∈ R and the diagonal matrices: Note, that (d 0 · p * ) corresponds to the global average density, (S -1 k · q * k ) to the local flow rate, and depending on the choices for p * and q * , the model's analytic and numerical character will differ.

Midpoint discretization
In case of the midpoint discretization, we set p * k and q * k to the mean of its associated edge's endpoints: Furthermore, we define q -:= 1 2 (q Rq L ), and note, that in vectorized form, p + = 1 2 (|A T 0 |p + |B T s |s p ). Together with the algebraic constraints from Sect. 2.3, a DAE system in the variables p, q + , and qarises: Since we aim to obtain an ODE, we need to transform this DAE system. The complexity of deriving this transformation is quantified by the DAE's index. From the various DAE index concepts, we use the tractability index τ [92], for which the midpoint discretization guarantees τ ≤ 2 [53].
This DAE can be decoupled into an ODE by rewriting it in the variables p and q + . To this end, 1. the pressure boundary condition implicitly resolves (2).
Since D p is a diagonal matrix, this is also numerically feasible.
We also pre-multiply (3b) with the inverse of the diagonal matrix D q . Altogether, we obtain: . (4b) This system of a pressure and mass-flux variable now consists of only differential equations. Notably, the first equation of the ODE system contains a temporal derivative of the input function s p , which practically would need to be approximated numerically, for example by finite differences. However, we will assume that all inputs are sums of step functions, so that effectivelyṡ p ≡ 0, which is reasonable as we assume exclusively lowfrequency boundary values.

Endpoint discretization
For the endpoint discretization, also called left-right discretization, we set p * k and q * k to the left and right endpoints, respectively: Since (B T s + |B T s |)s p = 0, we can write p R = A T 0,R p. With the algebraic constraints from Sect. 2.3, a DAE system in the variables p, q R , and q L results: As for the midpoint discretization, we want to derive a system of ODEs. For the endpoint discretization, it is shown in [52,115], that the tractability index is τ = 1, if all edges connecting supply nodes are directed from the supply, and each component of the graph is connected to at least one supply. This implies, no two supplies are to be directly connected. Similar to Sect. 2.4.1, this DAE can be decoupled into an ODE by rewriting it in the variables p and q L . Applying equivalent steps to (5a)-(5c) as for the midpoint decoupling (3a)-(3c) yields: using A 0,R + A 0,L = A 0 in (6a). An advantage of this endpoint discretization, in addition to the absence of a derivative of an input function, is the input-free friction term in (6b).

Temporal discretization
After spatial discretization and index reduction of the gas network model, a system of stiff nonlinear ODEs (in time) remains. The remaining temporal differential operator(s) are discretized using time-stepping schemes.
Due to the hyperbolicity of the pipeline model, the temporal resolution t used for the discrete integration interdepends on the spatial resolution x employed for the discretization of the spatial differential operators. Formally, this is expressed by the Courant-Friedrichs-Levy (CFL) condition [13,59,140], which states that the propagation of information cannot be faster than its conveyor: with the (dimensionless) CFL constant λ CFL , symbolizing the ratio of temporal and spatial discretization step-width scaled by the peak gas velocity v max . Since the flow is subsonic, v max could be estimated from the (linearized) characteristics [59], or via the boundary values. 6 However, we fix this maximum gas speed to v max = 20 m s (practically this is configurable in morgen).
Due to this relation of the space and time discretization and a pre-selected applicationspecific sampling frequency t of the output trajectory (i.e., every 60 s), x has to be adapted accordingly. The spatial discretization by finite differences in the previous sections ignores pipeline length, as each pipe is assigned only one (finite) difference. This means, pipes are potentially too long or too short with respect to a nominal length x, determined by the CFL condition x = (1ε)v max t, (0 < ε 1). Thus, too long pipes are subdivided into virtual pipes of nominal length, while too short pipes, including a potential remainder of too long pipes, are "rounded" to a full nominal-length pipe, yet with a friction term scaled by the fraction of the short pipe's length compared to the nominal length, This approach assumes that delays due to the forced virtual length of an actually short pipe are insignificant, hence this simple homogenization of pipe lengths may be improved by replacing short pipes with friction-less shortcuts and a static pressure-drop, as used in the quasi-static model [50,64] and similar to the subsequent compressor model in Sect. 2.5. Overall, we refine each pipe into a sequence of pipes of a selected nominal length -a graph level refinement -which is determined using the CFL condition. We note here, that this methodology is aimed at ensuring a certain minimum length for each pipe, as the shortest pipe may dictate an unnecessarily finely resolved time discretization. In terms of a pipe's maximum length, suggestions are for example: 5 km ( [50]), or 10 km ( [144]).
As the model is composed of a stiff, linear (hyperbolic) and a nonlinear component, an implicit solution of the linear part using a diagonally implicit Runge-Kutta (DIRK) method, and an explicit solution of the nonlinear part via a strong stability preserving (SSP) method, by an IMEX (IMplicit-EXplicit) solver, as proposed in [13,Sect. 3.2.2], is targeted. The actual quadrature rules used to compute the transient solutions are detailed in Sect. 5.3.

Simplified compressors
Beyond pipes, gas networks comprise a variety of non-pipe elements, of which the most important are compressors. Compressors increase gas pressure to counteract cumulative effects of retarding forces (friction, gravity, inertia etc.), and are grouped into stations with many possible configurations. For our purposes, we just allow fixed configurations on a macro scale [90] per scenario, which leads to a compressor being modeled as a special kind of edge that boosts the pressure from its suction inlet to its discharge outlet.
Compressors are typically modeled "ideally", based on their power consumption, for example as a special node type; in-depth discussions can be found in [63,123]. Due to this consumption model, such ideal compressor units are useful for energy utilization optimization tasks [40], yet, for a simplified transient simulation aspect already too complicated. A more practical approach is taken in [91,144,145], where a compression ratio α i ≥ 1 is used to scale the pressure in each node, or pipe [130], which means α i > 1 indicates compression / a compressor, otherwise (α i = 1) a pipe.
Here, we utilize a likewise simple compressor model similar to [46,139], for which we assume it is propelled by an external energy source, for example, given a compressor electrification, by excess renewable power [41], or that the off-take in gas is insignificant. But instead of using compression ratios, a constant (or parametric) target pressure is prescribed, modeling discharge pressure control [129].
The following affine compressor model is a variant of the compressor presented in [126] and used in [7]. We model a simplified compressor by a level, short pipe which increases the pressure at its outflow to a specified valuep c (and without friction, λ ij ≡ 0). Given the pipe from nodes i to j is treated as such a "compressor pipe", with target pressure ofp c , the variables p ij and q ij are given by the differential equations: The target pressurep c could be a step functionp c (t), and hence a control input [130], which could be accompanied by a discharge mass-flux output.
A compressor could also be interpreted as an actual pipe with "negative friction", and we considered using such nonphysical pipes as compressor model, but a difficult transformation between friction and pressure increase would have to be calibrated for every model variant (including friction factor formulas) and updated with every change in any model.

Parametrization
For the considered pipeline model (1), two scalar parameters are of practical interest: The temperature of the gas T 0 (in [K]), which is assumed to be constant throughout the network, and the global specific gas constant R S (in [J/(kg K)]).
Due to mainly underground on-shore pipelines [94,Ch. 45], and coolers in compressor stations [90], using an isothermal model is a reasonable simplification. However, temperature is relevant as a global parameter, since the use of an isothermal model "freezes" the dynamic energy (temperature) component in the original Euler equations (in time), and while intraday ambient temperature variation can be neglected for simulations with a 24hour time horizon, the temperature difference of a hot summer day and a cold winter day should be taken into account by a parameter representing an average temperature.
The specific gas constant, on the other hand, is determined by the gas composition, which may also vary. Again, local variations during an intraday simulation are neglected in this work, yet the average gas mixture of natural gas with, for example, hydrogen or biogas is relevant, so a parameter for the average specific gas constant is introduced. Together, the parameter-space is given by: and note that θ is used in the model (only) as Yet, lumping into a single parameter would impede physical interpretation.
Applied to the respective components of the input-output model in Sect. 2.7, this leads to parameter-dependent quantities, which need to be regarded accordingly by the model reduction as discussed in Sect. 3.4.

Input-output model
After spatial discretization and index reduction, we end up with a square input-output system, a system with the same number of inputs and outputs, consisting of an ordinary differential equation, an output function and an initial value: with parameter independent linear vector field components A and B, parametric mass matrix E(θ ), and nonlinear friction and gravity retarding term f (p, q, s p , θ ). The actual composition of the dynamical system components depends on the discretization and index reduction, cf. Sect. 2.4, while the linear output function C consists of C sq = -B s and C dp = B T d . The load vector F c ∈ R N q accumulates the respective discharge pressuresp c , as described in Sect. 2.5, for all compressors, and the initial state is given by a steady-state, whose computation is detailed in Sect. 2.8, depending on given steady-state boundary valuess p ,d q . Altogether, the gas network model is a generalized linear system (E, A, B, C), together with a nonlinear part f .

Figure 1 Schematic illustration of gas network input-output model
is justified by the numerical processing: only two model components depend on the parameters, temperature T 0 and specific gas constant R S , as well as on the compressibility factor z 0 , namely the mass matrix E p and the jointly treated retarding forces gravity and friction f q . Hence, the linear part of the right-hand side vector field is non-parametric and compressibility-agnostic. Overall, this system maps input boundary values, in the scope of this work, pressure at the inlets and mass-flux at the outlets, via the internal state, to output quantities of interest, here, mass-flux at the inlets and pressure at the outlets: To this type of input-output system we can now apply (data-driven) system-theoretic model reduction methods, which preserve the input-to-output mapping S, but explicitly not the internal state (p q) T . Lastly, we note that based on [74], we added a model fact sheet in the Appendix.

Steady-state computation
After spatial discretization, the dynamic simulation becomes an initial value problem. Yet, only the boundary values of the network model are known a-priori. This means the internal state at time t = 0 is unknown. We assume simulations always start at a steady-statē p,q for which ∂ t p = ∂ t q = 0, given some (initial) boundary valuess p ,d q . The internal state is then computable as a steady-state problem. Since the employed model is nonlinear, we approximate the steady-state by a two-step procedure: 1a. Linear mass-flux steady-state: 1b. Linear pressure steady-state: Step 2 can be repeated until an error threshold is met by using the previously approximated pressure steady-state. Practically, the linear problems in Step 1 and Step 2 are solved by a QR-based least-norm method [23]. Note, that Step 1a and Step 1b can be solved in parallel and that the QR decomposition of Step 1b can be recycled in Step 2 because of the chosen model structure. While this method works well for rooted-tree pipe-networks, it is not sufficient for cyclic networks with multiple supply nodes and non-pipe elements such as compressors. In this case, the resulting state after a limited number of the above algorithm's iterations is used as an initial value for the first order IMEX integrator detailed in Sect. 5.3.3, which time-steps until a steady-state is sufficiently approximated. This approximate steady-state, associated to a fixed set of boundary values and parameters, is used as initial value for the simulations: While other time steppers are applicable, too, the first order IMEX solver is related to the initial (two-step) algebraic approximation, due to the synthesis of the linear / input / source and nonlinear / reaction terms.

Port-Hamiltonian structure
An interesting class of models are port-Hamiltonian systems, which have already been used for gas network modeling [39,86]. Such port-controlled Hamiltonian models result from a system-theoretic approach to energy-based modeling, and are square, passive, stable and feature certain symmetries, besides their physical interpretability [15,103]. To exploit results from port-Hamiltonian theory in the context of data-driven model reduction, we regiment the previous modeling approach into the port-Hamiltonian framework. A linear input-state-output port-Hamiltonian model [135,Ch. 4] has the form: with a symmetric positive definite mass matrix E = E T , E > 0, a skew-symmetric energy flux J = -J T , a symmetric, positive, semi-definite energy dissipation R = R T , R ≥ 0, a symmetric, positive definite energy storage Q, Q > 0, resistive port matrix P and control port matrix G. 7 Here, we generalize the energy dissipation R ∈ R N×N to a nonlinear mapping R : R N → R N×N , this means the linear constraints become [38]: With this set up, we test the two index-reduced gas network model discretizations presented in Sect. 2.4.2 and Sect. 2.4.1 for compliance with the above port-Hamiltonian properties.
Proof We define the port-Hamiltonian state as x := (p q L ) T , which induces the remaining components. The mass matrix is symmetric positive definite, if its diagonal blocks are. Given that the D * are diagonal, and thus symmetric, as well as positive definite, both blocks are symmetric positive definite. The energy flux with the diag : R N → R N×N operator mapping a vector v to a diagonal matrix D such that v k → D kk , and element-wise (fraction) nonlinearities, results in one non-zero diagonal block and thus fulfills (10). The condition (11) is fulfilled since in the friction term of R, the absolute value of the mass-flux (numerator), and the nodal pressure variable (denominator) are always non-negative. Here, the energy storage represents the scale homogenization from Sect. 2.2, which is a diagonal matrix of positive entries, and due to same block structure in E also fulfills Q T E = E T Q. Lastly, the port matrix configuration complies to the port-Hamiltonian form.
Some remarks are in order on this result: From the previous proof it is also immediately clear that the midpoint discretization cannot be a port-Hamiltonian model, due to the input dependence of the energy dissipation. Furthermore, this derivation tests if the endpoint discretization has the mathematical port-Hamiltonian structure, but does not verify a physical energy-based model.
The somewhat nonphysical treatment of the gravity term as dissipating instead of storing ( [38]), is done with regard to the parametrization. Including the parametric gravity term as a retarding or damping force, and thus keeping the linear energy flux parameterfree, enables the previous steady-state computation.
Compressors, as modeled in Sect. 2.5 can be included by an additional summand inside the energy dissipation component, i.e. F C q , similar to the gravity term. This exhibits an unphysical negative sign inside the dissipation, as a compressor introduces energy. Furthermore, this compressor model requires to remove components from the A qp block of the system matrix (7), and thus perturbs the skew-symmetry of J.
Lastly, this notation for the dissipation can also be used for linearization, by constraining the argument of R to the steady-statex, Given the port-Hamiltonian model with a nonlinear resistive term, an (approximate) adjoint system can be derived by treating R as its image -a diagonal matrix. Transposing the (primal) port-Hamiltonian system's (9) transfer function , and exploiting the system components properties, yields the dual system: Hence, for the (nonlinear) endpoint discretization, its observability can be (approximately) quantified by the dual system's reachability, as for linear systems. Conceptually, this could also be done with the midpoint discretization, as it supplies the same model components.
However, it has no theoretical justification, as a dual system may not be accessible for (general) nonlinear systems.

Model reduction for gas networks
In this section, we summarize the principal approach behind all presented model reduction methods that are extended and tested in this work. The structure of the model laid out in Sect. 2 is given by (8). For large (expansive) networks, the differential equations in p and q become high dimensional, which impedes their solution [54,Sect. 7] and hence repeated simulations of scenarios. The aim of model reduction is to reduce the dimensionality of the differential equations, by computing subspaces of the phase space on which the trajectories evolve suitably similar (with regard to the quantities of interest). Furthermore, the reduced order model shall have the same form as the original model, and since two physical quantities are (bi-directionally) coupled in this system, the model reduction for interconnected systems [120] approach is used, yielding reduced operators for each sub-system: , centered around the steady-state (pq) T . This structure preserving model order reduction was already used in [18,53] in the context of model reduction for gas networks, while the centering has been used in [7] for gas network simulation and in [65] for nonlinear model order reduction. In the following, the general ansatz to obtain these reduced quantities (denoted by·) is summarized.

Projection-based model reduction
The reduced order model is computed by projecting the high-dimensional dynamics evolving in the (coupled) pressure and mass-flux phase spaces (of dimension N p and N q ) to low(er)-dimensional subspaces (of dimension n p and n q ), which capture the principal components of the respective trajectories. Given suitable discrete projection mappings from the original space to the reduced space V T * and mappings from the reduced space back to the original space U * : Thus, the reduced trajectory results from applying V * to the original trajectory's steadystate deviation, while the original trajectory is approximately recovered by applying U * to the reduced trajectory: the initial condition is also reduced by application of V * . Similarly, the components of the reduced system result from applying the U * map to the argument of the respective operators, and the V * map to the result of the operation. For the linear operators, the matrices E * , A * , B * and C * and the vector F c , this leads conveniently to pre-computable reduced matrices and vector respectively, C dp := C dp · U p ∈ R N d ×n p , yet, the nonlinear componentf q remains a composition operation:

Structure preserving model order reduction
In this specific context, the term structure preserving model order reduction (SPMOR) has two meanings: first and foremost, SPMOR refers to the separate reduction of the state components, as above in the case of gas networks, the individual reduction of the discretized pressure p and mass-flux q variables. Second, SPMOR can also refer to preserving the port-Hamiltonian form (9). For projection-based model reduction, the former is generally ensured by separate projectors [44] (or an overall block diagonal projection). The latter is guaranteed by using Galerkin projections, which implies stability preservation [16], given a port-Hamiltonian full order model. Both SPMOR interpretations are jointly fulfilled if a block-diagonal (w.r.t. p and q) Galerkin projection is used.

The lifting bottleneck and hyper-reduction
The gas network models considered for reduction are nonlinear (and potentially nonsmooth), hence the reduced order nonlinear partf q involves lifting the reduced state up to the original high-dimensional space, evaluating the nonlinearity and projecting the result back down to the reduced low-dimensional space (13). As the high-dimensional space is involved, this is typically computationally demanding and may eat up the gains from the reduction of the linear part. To mitigate this so-called lifting bottleneck, hyper-reduction methods can be employed, which construct low-dimensional surrogates for nonlinearities.
In this work we discard (or rather defer) hyper-reduction due to the following reasoning: The purpose of this work is to determine which method constructs the best reduced order models, a hyper-reduction may inhibit comparability due to, for example, a dominating hyper-reduction approximation error. Second, various hyper-reduction methods for this setting are applicable (i.e. DEIM [26], Q-DEIM [36], DMD [142] or numerical linearization [98]), which may interact differently with the different model reduction methods. So as a first step, the bare model reduction methods are tested here (this means: which method's linear subspaces capture the nonlinear dynamics best), at a later stage the best hyper-reduction method can then be determined. Lastly we note, the nonlinear part of the vector field consists exclusively of element-wise operations (see Sect. 2.9), a system with repeated scalar nonlinearities (SRSN) [28], which are less difficult to handle due to "locality" of the nonlinearity, and hence, its vectorization.

Parametric model reduction
There are two common approaches for projection-based parametric model order reduction: averaging and accumulating [67]. For the selected data-driven methods, averaging means that for a set of parameter samples the associated trajectories or derived quantities (such as the utilized system Gramians) are averaged, while accumulating refers to the concatenation of trajectories or derived quantities (such as the projectors).
Generally, each of the structure-preserving model order reduction methods in Sect. 4 can be used with either, we opted to use the averaging ansatz for all of the following methods since their computation is without exception based on parametric empirical Gramians [71].

Model reduction methods
In this section, we briefly summarize the employed model reduction methods from a practical point of view. For theoretical details and backgrounds we refer to the relevant works, cited in the respective subsections. Due to the non-differentiable nonlinearity (friction), the sought projections U * , V * for all tested model reduction techniques are constructed from (transformed) time-domain trajectory data obtained from numerical simulations, which is given by discrete-time snapshots of the internal pressure nodes X p and mass-flux edges X q , the external demand node pressure Y p and supply node mass-flux Y q , , as well as dual state components Z p and Z q , if available, The state-space trajectories X * (t; θ k ), Z * (t; θ k ) are obtained for perturbations of the inputs, pressure at the N s supply boundary nodes and mass-flux at the N d demand boundary nodes, while the output trajectories Y * (t; θ k ) are computed for perturbations in the respective N * steady-state components. Using the dual state trajectories is significantly faster than output trajectories, as computing observability as dual reachability scales, as for the primal reachability, with the number of ports, instead of scaling with the number of internal states. The training parameters θ k are sampled from a sparse grid spanning the parameter space ⊂ R 2 . All methods are prefixed "Structured", since the model structure of a pressure and massflux variable is preserved in the reduced order model. Practically, this means while pressure and mass-flux trajectories are computed simultaneously due to their coupling, the individual projectors for pressure and mass-flux are constructed separately.
The subsequent methods may not have been previously introduced explicitly in structured form, yet given [1,120,136] introducing structured Gramians, these are trivial extensions. For ease of notation, we describe the computation of projectors {U p , V p } and {U q , V q } generically as {U * , V * }.
We implemented a total of thirteen model reduction method variants, which we compare in this work. All tested methods are data-driven and time-domain focused, as the dynamic gas network model (8) is nonlinear. Furthermore, all model reduction methods construct linear subspaces and are derived from methods for linear systems, yet differ from plain linearization: instead, the implemented methods assemble linear subspaces of the model's phase space that, in a method specific sense, approximately enclosing the relevant nonlinear system evolution. Moreover, all methods are SVD-based [4], and their majority is originally based on (empirical) system Gramian matrices, for details see [66].
We highlight here that the time horizon for the training data is significantly shorter than for the actual simulations the reduced order model is targeted at. Furthermore, generic training inputs (transient boundary values), such as impulse, step or random signals, are utilized to avoid a model reduction crime [67] (comparable to an inverse crime): Test a reduced order model using the training parameters or inputs.

Empirical system Gramians
All model reduction methods currently included in morgen are computationally realized using empirical system Gramian matrices, which are system-theoretic operators encoding reachability and observability. From these, information on importance of linear combinations of states can be extracted. For linear systems, these system Gramians are typically computed via matrix equations, for general nonlinear systems there is no feasible closed form. However, the empirical system Gramians approximate the nonlinear Gramians based on state and output trajectory data. In case of a port-Hamiltonian (nonlinear) systems, the approximate dual system (12) enables substituting expensive state by port perturbations, and thus severely reduce computation times. Empirical Gramians are described in-depth in [66]; following only a brief summary is given.

Empirical reachability Gramian matrix
Reachability quantifies how well a system can be driven by the inputs, which is encoded by the reachability Gramian. The empirical reachability Gramian is an approximation based on state trajectory data [84]: Given a suitable set of input perturbations, W R, * approximates the nonlinear reachability Gramian near a steady-state [56].

Empirical observability Gramian matrix
Observability quantifies how well the state can be characterized from the outputs, which is encoded by the observability Gramian. The empirical observability Gramian is an approximation based on output trajectory data [84]: For a suitable set of steady-state perturbations, W O, * approximates the nonlinear observability Gramian near a steady-state [56]. Given the port-Hamiltonian structure of a discretization, the empirical reachability Gramian of the dual system (12) can be used to compute the empirical observability Gramian [67,Sect. 2] with Z m * (t; θ k ) instead of X m * (t; θ k ).

Empirical cross Gramian matrix
The empirical cross Gramian concurrently encodes reachability and observability, which in conjunction quantifies redundancy also known as minimality; however, the (empirical) cross Gramian is only applicable for square systems, systems with the same number of inputs and outputs (which the gas network model (8) fortunately is), and an approximation is based on simulated state and output trajectory data [70]: For a suitable set of input and steady-state perturbations, W X, * approximates the nonlinear cross Gramian near a steady-state [70]. Given the port-Hamiltonian structure of a discretization, the linear empirical cross Gramian exploiting the dual system (12) can be used to compute the empirical cross Gramian [14] with Z m * (t; θ k ) T instead of Y * (t; θ k ).

Empirical non-symmetric cross Gramian matrix
A generalization of the empirical cross Gramian for non-square systems is the empirical non-symmetric cross Gramian, which is an approximation based on simulated state and (averaged) output trajectory data [72]: For a suitable set of input and steady-state perturbations, W Z, * approximates the nonlinear cross Gramian near a steady-state. Even though the gas network model (8) is square, the empirical non-symmetric cross Gramian is included here, since, heuristically, it could provide better results than the regular cross Gramian [72]. Furthermore, an empirical nonsymmetric linear cross Gramian is computable by similarly averaging over the dual states and replacing Y q * (t; θ k ) by Z q * (t; θ k ) T .

Structured proper orthogonal decomposition
Proper orthogonal decomposition (POD) is a basic data-driven method for model reduction: given a matrix of state snapshots over time, the dominant left singular vectors are computed as a basis via a singular value decomposition (SVD) [97]. The basis vectors are assigned their principality with respect to the conveyed energy by the associated (relative) singular value magnitude.
In the context of this work, the POD is constructed from a system-theoretic point of view, that connects to the system property of reachability. Due to the overall structured approach to model reduction, a structured POD refers in this context to the separate PODs for pressure and mass-flux variables p and q as in [18,51,53].

Reachability-Gramian-based
The singular vectors to the principal singular values of the empirical reachability Gramian (14) correspond to the POD modes. To obtain a reduced order model, first, a truncated SVD (tSVD) of the empirical reachability Gramian, reveals the principal subspace of the respective trajectories X m * (t), whereas the importance of each basis vector (column) in U R, * is determined by (the square-root of ) the associated singular value σ i := D 1/2 R, * ii : The matrix of basis vectors (POD modes) constitutes a Galerkin projection V * := U * . Notably, (structured) POD only considers the input-to-state mapping, not the state-to-output mapping, and hence approximates the state variables, p and q, not the output quantities of interest s q and d p .
The POD could also be computed directly from an SVD of the trajectory data, yet the computational overhead of using the empirical reachability Gramian is small compared to the trajectory simulation runtimes, and the systematic perturbations of the empirical Gramian approach [66] are exploited.

Structured empirical dominant subspaces
The previous (structured) POD method considers only the reachability information, hence the data only reflects the input-to-state mappings, and thus the POD derived ROMs (Reduced Order Model) approximate the state variables p and q. To approximate the outputs s q and d p , the state-to-output mappings, encoding observability information, need to be considered, too.
The (empirical) dominant subspaces method initially developed in [110], and originally named DSPMR (Dominant Subspace Projection Model Reduction), conjoins and compresses the dominant reachability and observability subspaces of an input-output system, such as the gas network model (8), obtained from (empirical) system Gramians. Heuristically, this method seems to be useful for hyperbolic input-output systems [49]. Here, we consider three variants: first, based on the empirical reachability and observability Gramians, second, based on the empirical cross Gramian and third, based on the empirical nonsymmetric cross Gramian.

Reachability-and observability-Gramian-based
The singular vectors associated to the principal singular values of the empirical reachability and observability Gramians span these dominant subspaces, which are first extracted by tSVDs, and then, after concatenation, compressed by orthogonalization via another tSVD: F equilibrate the potentially mismatched scales of the respective Gramians, akin to the refined DSPMR method from [110].

Structured empirical balanced truncation
The dominant subspaces method combines separately quantified input-to-state and stateto-output energies, but not the actual input-to-output energy. Such can be accomplished by balanced truncation and based on the Hankel operator, which maps past inputs to future outputs. This operator's singular values measure the sought input-to-output energy, and the singular vectors constitute a basis. To obtain the Hankel operator's truncated SVD, first, the underlying system needs to be balanced, and then singular vectors associated to small magnitude Hankel singular values are truncated.
Balanced truncation is the reference model reduction method for linear input-output systems, due to stability preservation in the ROMs and a computable error bound. For (control-affine) nonlinear input-output systems, such as the gas network model (8), balanced truncation can be generalized to empirical balanced truncation [84], while a structured variant is introduced as interconnected system balanced truncation [136], which we combine.
Again, we consider three variants: First, based on the empirical reachability and observability Gramian, second, based on the empirical cross Gramian and third, on the empirical non-symmetric cross Gramian; additionally, the reachability and observability-based balanced POD variant is included.

Reachability-and observability-Gramian-based
The original balanced truncation method is based on the reachability and observability Gramians [99]. A transformation into a balanced coordinate system in which both system Gramians are diagonal and equal is obtained via simultaneous diagonalization. Various balancing algorithms are available for this task [137], in this setting we selected the general balancing algorithm [118,119], which utilizes the magnitude-based truncated eigenvaluedecomposition (tEVD) of the Gramians: The matrices U * , V * constitute a Petrov-Galerkin

Cross-Gramian-based
For linear, symmetric systems, alternatively a balanced and truncated reduced order model can be computed via the cross Gramian. Yet, the gas network model (8) is neither linear nor symmetric, not even in a nonlinear sense of symmetric systems, i.e. gradient systems [75], but the considered input-output system is square. Hence, an empirical cross-Gramianbased reduced model is computable, but it will differ from the reduced model obtained by reachability-and observability-Gramian-based balanced truncation. Given a cross Gramian with full rank, an (approximate) balancing projection is computable in a similar manner as for the (empirical) balanced truncation (18), but based on the left and right eigen spaces of the cross Gramian [77]: The matrices U * , V * again constitute a Petrov-Galerkin projection, and the importance of each column is determined by the absolute value of the (diagonal) elements of D X, * , which are only equal to the Hankel singular values for linear and symmetric systems. As for the empirical dominant subspaces method, the cross-Gramian-based empirical balanced truncation variant directly extends to the non-symmetric cross Gramian W Z, * .
One could assume that if only one system Gramian has to be computed, instead of two for dominant subspaces or balanced truncation, that the cross-Gramian-based computation is significantly faster, but the overall number of simulated trajectories is the same for both methods, which causes the dominant fraction of computational cost. Thus, the empirical cross Gramian computation is merely somewhat quicker than the computation of both empirical reachability and observability Gramians.

Structured balanced POD
Instead of balancing the system using the product of the (empirical) reachability and observability Gramians, the dominant subspaces of the respective Gramians can be used to approximately balance the system. This variant of reachability-and observability-Gramian-based balanced truncation is called balanced POD [141], and the approximate balancing algorithm reads: with the matrices U * , V * inducing a Petrov-Galerkin projection.
Here, we categorized balanced POD as a variant of balanced truncation, due the algorithmic similarity to the balancing algorithm (18). Alternatively, we could have classified balanced POD as a variant of POD, since it can be described as POD with the observability Gramian defining the POD's inner product [116].
We leave it to the reader to test other balancing algorithms, e.g. [17,137], while we excluded the modified POD [67,102], as it is not a Petrov-Galerkin method.

Structured empirical balanced gains
The empirical balanced gains method is a variant of the empirical balanced truncation method [30,67]: While balanced truncation selects principality of subspaces based on the Hankel singular value magnitude, balanced gains [30] sorts the balanced basis vectors by the impulse response norm: 2 ≈ tr C dp W R,p C T dp ≈ tr B T pd W O,p B pd ≈ tr(C dp W X,p B pd ) , for a system in balanced form. Hence, either method, balanced truncation and balanced gains compute the same balancing transformation, via (18) or (19), but the sequence of basis vectors differs due to the variant measures. We note, that due to the linear input and output operators, the balanced gains method can be directly applied for the nonlinear gas network models (8).

Structured goal-oriented proper orthogonal decomposition
Similar to the balanced gains method, (structured) POD basis vectors can also be sorted, instead of by their singular value magnitude σ k , in terms of their impulse response, akin to the simplified balanced gains in [30], by the index d k for rowsc k of the POD-transformed output matrix C. Given the reachability-based POD, this variant is related to the concept of output-reachability and H 2 -norm model reduction [112].

Structured dynamic mode decomposition -Galerkin
In addition to the four energy-based method classes (Sects. 4.2, 4.3, 4.4, 4.5), also an alternative method from [2], based on system identification, is investigated. Yet, it is adapted as a structured projection-based model reduction method. After summarizing the parent system identification method, the derived model reduction technique is presented.

Dynamic mode decomposition
Dynamic mode decomposition (DMD) identifies a discrete-time operator from (discrete) trajectory time-series data, preserving certain modal behavior [117]. This system identification method is based on the Koopman operator, which is an infinite dimensional, but linear operator, mapping (a transformation, or observable of the) state x k at discrete-time step k to (k + 1). DMD yields a (linear) finite-dimensional approximation of the Koopman operator preserving its dominant eigenmodes. Here, using the identity observable and given time series data X = [x 0 . . . x K ], DMD identifies an operator A, via least-squares: Since the underlying model is a control system, DMD with control (DMDc) [113] (and known input operator B), i.e. using X c := [x 0 -Bu 0 . . . x K -Bu K ] instead of X for DMD, is applicable, yet, due to the perturbed steady-state training of ROMs (see Sect. 6.1) not beneficial.

Reachability-Gramian-based DMD-Galerkin
DMD is rather a system identification than a model reduction method. To fit into the projection setting, we utilize the DMD-Galerkin method [2], which forms a Galerkin projection U from the orthogonalized dominant eigenvectors (based on the associated eigenvalue magnitudes), Since only state-space trajectories are utilized, as for the POD, the DMD-Galerkin method approximates the state, not the output. Curiously, we note, that this method uses discretetime information to assemble projections for continuous-time systems. Practically, U * are computed as a subset of singular vectors for the largest magnitude singular values, instead of orthogonalized eigenvectors. In the structured setting at hand , with X * representing the pressure and mass-flux states X p , X q , the respective Galerkin projection U * = V * is given by: Effectively, the approximate Koopman operator A is computed via an empirical reachability Gramian, yet instead of the standard inner product, the DMD-"kernel" (21b) is used. Additionally, and in line with the original empirical Gramians [84], the utilized trajectories are centered, following [73]. This means theoretically, DMD-Galerkin is POD with a specific kernel, and practically, that by computation via an empirical Gramian the systematic perturbation properties can also be exploited for DMD-Galerkin.

morgen
The morgen (Model Order Reduction for Gas and Energy Networks) platform (version 1.0) implements the mathematical methods presented above, in MATLAB (≥ 2020b) and compatible to Octave (≥ 6.1). Compared to, for example [62], morgen does not feature a graphical user interface or Simulink integration, since it is designed for batch testing and multi-query use on (headless) workstations.
The source code is organized into five main components: networks (holds network and scenario data-sets) models (discretizes networks and assembles input-output systems) solvers (computes solution time series for discrete models and scenarios) reductors (reduces state variables of discrete models) tests (defines experiments for data-model-solver-reductor combinations). These components are briefly described in the following. For an illustration of the internal structure of morgen, see Fig. 2. Further code that is used by the main function or multiple components is contained in an utils folder, and the stand-alone network format converters 8 are stored in a tools folder.

Networks
The networks directory stores the network .net files, and a folder for each network with the same name as the associated network file's base name, which hold the scenario definition files. The .net file, in comma-separated value (CSV) format, defines the network topology and gas network components through an edge list. Each row encodes one edge through the information: edge type (pipe, shortcut, compressor, valve), "from"-node, "to"-node, length, diameter, incline, roughness; the latter three are only relevant for pipe edges.
Note that boundary nodes have to be leaf nodes of the network graph in order to be identified by morgen as such. Furthermore, if a boundary node shall act as supply and Figure 2 Internal data flow and process of morgen demand, it needs to be artificially split into two leaf nodes, as the edge connecting a supply node has to be directed from the leaf, while the edge connecting a demand node has to be directed towards the leaf (Sect. 2.3).

Scenario
A scenario is defined via a set of key-value pairs in an .ini file, whereas the sequence of pairs does not matter. Each network has at least a generic training.ini scenario. The following keys are mandatory for a scenario definition: T0 Mean temperature RS Mean specific gas constant tH Time horizon ut List of time instances up List of supply node pressures at time instances ut uq List of demand node mass-fluxes at time instances ut. Depending on the network composition, the following keys may need to be provided: cp List of compressor discharge pressures vs List of valve settings. Note that amongst other configurations, the parameter ranges of temperature and specific gas constant for the parametric model reduction are set in the global morgen.ini file.

Models
A model encodes a spatial discretization of the simplified Euler equations (8) on a gas network topology in a structure with the members: A, B, C, E, F, f and the Jacobian J ( [3,128]), together with the system dimensions in terms of number of pressure-at-nodes and mass-flux-on-edges states, and total number of boundary / port nodes. It is ensured during the assembly of the model that the sparsity of the model components is preserved.
The model interface is given by the following signature: discrete = model(network,config); While the linear components A, B, C, and F are provided as sparse matrices, the parameter-dependent linear component E is a closure 9 returning a sparse matrix, and the nonlinear component f and the (nonlinear) Jacobian J are closures returning the application of a state (as well as steady-state, input, parameters, compressibility).
Two spatial discretizations are currently provided in morgen: Even though both provided models are ODEs, DAE models can also be implemented and tested, given a solver (and reductor) is available.

Solvers
The morgen platform provides four solvers: An adaptive step-size method, a fixed stepsize explicit method and two fixed step-size implicit-explicit methods.
While explicit solvers only require vector field evaluations at the cost of smaller time steps, implicit solvers have to solve a root-finding problem in each time step for a nonlinear model. An IMEX with singly DIRK (SDIRK) solver turns out to be the most efficient for this class of models, simplifying the nonlinear problem to a linear problem solvable by a single matrix decomposition per trajectory. Due to the non-diagonal mass matrix, even in the case of a sufficiently stable and accurate explicit method, at least one linear problem per trajectory would have to be solved.
The solver interface is given by the following signature: solution = solver(discrete,scenario,config); with the return value solution, being a structure, and the arguments discrete (model), scenario, and config(uration).
All provided fixed step-size solvers cache matrix decompositions. The initial steadystate is also cached, as is the QR decomposition used to compute it.
Even though the overall model (8) has a two dimensional structure and the reductors exploit this structure, the simulations itself can be performed on a lumped model (we omit parametrization here for ease of notation):

Second-order adaptive implicit solver: generic
For validation purposes, the adaptive step-size solver for stiff systems ode23s 10 included in MATLAB (and Octave) [125,Sect. 2.3], based on a modified second order Rosenbrock formula is used and encapsulated as a generic solver. Due to preferential performance demonstrated in [51], ode23s is preferred over alternatives such as ode15s (or ode45).

Fourth-order "classic" Runge-Kutta solver: rk4
Since in [104,105], the fourth-order explicit Runge-Kutta (RK) method [83] is employed, morgen provides it, too. This method is strong stability preserving [47], however it is not SSP-optimal and it works only for small time-steps.

First-order implicit-explicit solver: imex1
The lumped gas network model (22) can be split into a linear and a nonlinear part, of which the linear part is numerically stiff and hence should be solved with an implicit solver, while an explicit solver is preferred for the nonlinear part as to avoid solving a root-finding (optimization) problem in each time step.
An IMEX method allows this separate treatment of operators and is thus suitable for this hyperbolic and nonlinear system. Combining the first-order explicit Euler's method with the first order implicit Euler's method yields the first order IMEX method: with the associated Butcher tableaus Table 2. Even though this IMEX method is not a Runge-Kutta method [5], it was successfully applied to gas network models in [52] with a relaxation parameter set to γ = 1.

Second-order implicit-explicit Runge-Kutta solver: imex2
A second-order (two-stage) IMEX Runge-Kutta solver is provided, based on the combination of a second-order explicit SSP Runge-Kutta method [47], and a second-order DIRK method. Following [109], such an IMEX-SSP2(2,2,2) method with relaxation γ is given by: The explicit component of this IMEX method is SSP-optimal [47], while depending on the choice for the free parameter λ, different properties of the implicit component can be achieved (see Table 3). Practically, we found λ = 1 2 , making the implicit part SDIRK and stiffly accurate, to work best. Additionally, we would like to highlight passive Runge-Kutta methods [101], specifically PDIRK (passive DIRK), as implicit IMEX component, which have various desirable stability and conservation properties [43].
The associated Butcher tableaus are given in Table 4.
In our experiments, the first-order IMEX method Sect. 5.3.3 allowed larger time-steps and exhibited less numerical oscillations or artifacts compared to the second-order IMEX-RK methods. The generic (adaptive) method Sect. 5.3.1 also solves sufficiently accurate, but takes longer to compute. Thus, by default, we recommend the first-order IMEX integrator for gas network simulations.

Reductors
The reductor module provides methods that compute structured 11 (struct.) projectors for a given discretization. These projectors can be stored on disk for reuse. Currently, the reductors, described in Sect. 4 are included: Each reductor variant has a suffix characterizing the employed empirical Gramians: _r for only reachability, _ro for reachability and observability, _wx for cross, and _wz for non-symmetric cross Gramian. For each reductor utilizing observability information (this includes the cross Gramians), a linear variant using the dual system is available, and signified with the additional suffix _l.
The reductors have the interface: [proj,name] = reductor(solver,discrete,scenario,config); returning a (cell) array of projectors with maximum configured column rank, as well as the reductor's full name. These projectors specific to model and solver defining a ROM are then stored in a .rom file.

Empirical Gramian framework
The compute back-end for the model reduction methods is emgr -empirical Gramian framework [66]; currently in version 5.9 [68]. This (open-source) Octave and MATLAB toolbox computes the empirical Gramians, which are essential to construct the reduced order models via the Gramian-based model reduction methods from Sect. 4, including the structured DMD-Galerkin method.

Tests
A test is a script defining an experiment by specifying network, scenario, model, solver and reductors. The tests component is a collection of test scripts probing primarily model reduction for various networks. A typical test contains two calls to the main morgen function. The first call computes the reduced order model (offline phase): morgen(network,training_scenario,model,solver,reductors); which computes a ROM from a short, generic, steady-state training_scenario. The projectors defining the reduced order models are then stored in rom_files. The second call tests the reduced order model(s) on a longer test_scenario (online phase): morgen(network,test_scenario,model,solver,rom_files); Generally, a model reduction method can also be tested in a single call, disregarding that scenario's boundary value time series, yet for productive use, a reduced order model is constructed once (first call), and then employed for many different scenarios (second call).
Included in morgen are two types of tests: First, tests prefixed with "sim_" only simulate the test scenarios, second, tests prefixed with "mor_" compute the reduced order models using training scenarios, and benchmark these ROMs on the test scenarios.

Numerical experiments
In the following, we present three sets of numerical experiments, with the purpose of demonstrating the reducibility of gas network models via the data-driven, parametric, system-theoretic model order reduction algorithms from Sect. 4, and illustrating the capabilities of the morgen framework summarized in Sect. 5. The first set uses a pipeline "network", which is interesting in the context of model reduction, while the second set tests an academic toy network as a sanity check and a simple functionality test. Lastly, a realistic gas network topology is evaluated.
We note that various further networks are included for testing in morgen; among others: the Canvey-Leeds network [55,79], the Belgium transport network [31,95] and a part of the Fermaca network [114]. A synthetic pipeline model and associated simulation results were provided by the PSI Software AG for validation of morgen against the commercial PSIganesi 12 solver.

Workflow
For each of the numerical tests, the same workflow is employed, which is composed of a training phase (offline phase), in that the ROMs are computed, using a generic test scenario, with a (virtual) time horizon of 1 h, and boundary value input functions typical for system identification, i.e. Dirac impulse, step signal, random-binary signal or Gaussian noise [100,Ch. 16]. In the test phase (online phase), the ROMs are tested on scenarios with a (virtual) time horizon of 24 h [27,90], (starting at 6 am [40]). In addition to shorter offline phases, this difference in training and test time horizons emphasizes generality of the ROMs.
To verify models and solvers [107] this offline/online procedure is performed for all combinations of: whereas the port-Hamiltonian ode_end model is tested with the nonlinear as well as the linear reductor variant (_l suffix) if available, while the ode_mid model is only tested with the nonlinear reductor variant. The models are specialized by the Schifrinson friction, and the AGA88 compressibility factor formula. We excluded the generic and rk4 solvers in this comparison as they are too slow or too fragile, respectively. Yet, the test scenario visualizations in Fig. 4(a) and Fig. 6(a) are computed by the generic solver.
For the parametric model reduction, the temperature range for training and testing is set to [0 • , 15 • ]C, while the specific gas constant range is chosen as [500, 600] J kg K . During training, samples from the parameter space are drawn from a sparse grid, whereas for the tests, parameters are drawn from a uniform random distribution. For either test and training, five parameters are sampled. The input perturbations for the steady-state training scenario are selected to be a step function, which heuristically works well for hyperbolic systems [49]. The reduced order models are compared via the approximate, discrete, parametric, (L 2 ⊗ L 2 )-norm of the output error [67]: for a finite sample h of the parameter space and discrete output samples y h (θ h ),ỹ h (θ h ). This energy norm is chosen, since all methods are at least related to an energy-based method. However, morgen can also provide the errors in the approximate parametric (L k ⊗ L ) parameter-space-state-space norms for k ∈ {1, 2, ∞}, ∈ {0, 1, 2, ∞}, cf. [50]. Note, that due to nonlinearity of the considered models, and the averaging nature of the norm, a monotonic error decay cannot be expected. To enhance comparability of the results, also the reducibility measure MORscore [67] for each experiment is computed. The MORscore for a certain method and model is essentially the area above the method's error graph in the relative error plot such as Fig. 4(c).
Lastly, we note that the following numerical experiments are conducted using a computer with an AMD Ryzen 4500U @ 2.3Ghz hexa-core processor and 16 GiB memory running MATLAB 2021a on Ubuntu 20.04 Linux.

Yamal-Europe pipeline
First, a pipeline is tested, which is an interesting test case, since the trivial topology (Fig. 3) comprises little redundancy, hence pipelines are a useful benchmark for model reduction methods.
The Yamal-Europe pipeline connects gas fields in Russia with western Europe. 13 A section of this pipeline was also benchmarked in [19,24,25,107], from which the technical properties and test scenario are taken. The considered pipeline section is 363 km long, has a diameter of 1.422 m, no (reported) inclination, and a pipe roughness of 0.01 mm. A steady-state, used as initial state, is set by a supply pressure of 84 bar and demand massflux of 46.3 kg s . The semi-discrete nonlinear state-space system has two inputs and outputs as well as 908 states; and a time step width of 20 s is used. The employed test scenario is taken from [25], compressed to 24 h, and shown in Fig. 4(a), the associated model reduction errors are given in Fig. 4(c), Fig. 4(d), Fig. 4(e), and Fig. 4(f ), for up to reduced order 150, while the resulting MORscores are listed in Table 5.
Generally, the choice of solver is more relevant than the choice of model: while the MORscores for different models but same solver are similar, for the same model but different solver, they are significantly dissimilar. Also, the tested balancing (Petrov-Galerkin) methods perform worse than the Galerkin methods.  For both models, and the first-order IMEX solver, the structured empirical dominant subspaces methods perform best, followed (closely) by the DMD-Galerkin and (goaloriented) POD method; then, among the balancing methods, the balanced POD and cross-Gramian-based variants. The most overall accurate method is the cross-Gramian-based dominant subspaces method.
For both models, in combination with the second-order IMEX-RK solver, the structured POD and goal-oriented methods lead, followed by the DMD-Galerkin reductor. For both solvers, the endpoint model performs better than the midpoint model. In case of the port- Table 5 MORSCOREs μ(150, mach (16) ) in the L 2 ⊗ L 2 error norm for the "Yamal-Europe" pipeline network from Sect. 6 We note that the cross-Gramian-based dominant subspaces methods produce the lowest errors, and since for the linear reductors used in combination with the endpoint model, the dominant subspaces methods are as efficient as the purely reachability-based DMD-Galerkin, and (goal-oriented) POD methods.
Interestingly, while the second-order IMEX-RK solver is better suited for simulations of the full order model, in terms of data-driven model reduction and/or reduced order model simulation it is significantly worse than the first-order IMEX solver. This is also demonstrated in the subsequent experiments.

MORGEN network
The second set of tests encompasses a synthetic network for testing morgen's capabilities. This "MORGEN" network, with topology as in Fig. 5, tests the interaction of simplified compressors from Sect. 2.5 for various network features, such as cycles, multiple supply and demand nodes, and is in the spirit of a test network from [39].
Specifically, six sub-networks (in the shape of letters) are connected, the second and third sub-network contain a cycle, a compressor connects the third and fourth subnetwork, and the fourth and fifth sub-network contain additional supply and demand nodes. The edges vary in length between 20 km and 60 km, while the diameter and roughness are consistently 1 m and 0.01 mm, respectively. A steady-state, used as initial state, is set by supply (and discharge) pressures of 50 bar at both supply nodes and the compressor, and demand mass-fluxes of 30 kg s at all demand nodes. In semi-discrete form, the nonlinear state-space system features six inputs and outputs as well as 901 states; and a time discretization with 60 s time steps. The employed 24 h test scenario is made from hourly standard load profiles [60] and shown in Fig. 6(a) Table 6. Again for this comparison, the choice of solver is more relevant than the choice of model, yet the MORscores are much lower, due to the complexities (cycles, compressor, multiple demands) of the network.
Both models in conjunction with the first-order IMEX solver only produce workable results with Galerkin methods. The endpoint ROMs are again more accurate than the midpoint ROMs. Notably, the linear reachability-and-observability dominant subspaces method for the endpoint model, is leading the MORscores.
The ROMs for both models with the second-order IMEX-RK solver perform worse, with the exception of the balanced truncation ebt_ro and balanced gains ebg_ro variants. However, the second-order IMEX-RK solver related ROMs are of no practical use due to the high(er) error.
Overall, the (linear) eds_ro dominant subspaces reductor produces the lowest error, followed by the DMD-Galerkin, (goal-oriented) POD, and cross-Gramian-based dominant subspaces methods. As for the pipeline, the endpoint model is better suited for model reduction, while the first-order IMEX solver results in significantly more accurate ROMs than the second-order IMEX-RK solver.

GasLib network
Lastly, a network topology derived from real-life is tested. The GasLib-134v2 network [122], modeling a part of the Greek natural gas transport system, is overlayed on a map of Greece in Fig. 7. The network has a total length of 1412 km and features a compressor. A steady-state, used as initial state, is set by supply (and discharge) pressures of 80 bar at supply nodes and the compressor, and demand mass-fluxes up to 16 kg s at all demand nodes.
In semi-discrete form, the nonlinear state-space system has 48 inputs and outputs as well as 2682 states; and 30 s time steps are employed. For testing, a random (24 h) load profile is generated, by adding samples from a scaled uniform random distribution to the steady-state, 14 shown in Fig. 8(a), the associated model reduction errors are given in  Table 7.
As before, the choice of solver is more relevant than the choice of model. Challenges in this network, beyond the compressor, are the high number of boundary nodes, which are predominantly demand nodes (N d = 45).
First, we note that only Galerkin methods produce consistently stable ROMs. Furthermore, in comparison with the previous experiments, the dominant subspace methods perform worse, and all variants based on reachability and observability Gramians perform relatively better. The endpoint model seems to be better suited for the tested model reduction methods than the midpoint model. And as for the other experiments, the first order IMEX solver outmatches the second order IMEX-RK solver.
Considering all experiments, the DMD-Galerkin method performs best in terms of MORscore, accuracy and efficiency, followed by the dominant subspaces methods. We also note that the purely reachability-based as well as the linear reductors exploiting the port-Hamiltonian structure are the most efficient. Surprisingly, Galerkin methods perform better than the tested Petrov-Galerkin methods in terms of accuracy and stability, while in an unstructured, non-parametric, linearized setting all tested Petrov-Galerkin methods would be stability preserving. Yet, structured balancing methods are explicitly not guaranteed to be stability-preserving [120,136].
With regard to the computational complexity of the offline and online phase, we reiterate, that due to the absence of hyper-reduction, the online runtimes are not competitive (see Sect. 3.3), thus, we focus on the offline phase. Yet, due to the practical reducibility of the state-space dimension by more than one order of magnitude in the numerical experiments using the first order IMEX solver, a considerable speed-up is to be expected.
For the tested data-driven (time-domain) model reduction methods, the number of vector-field evaluations, or relatedly, the number of simulated trajectories measures the complexity, as these constitute their principal fraction. The empirical reachability Gramian requires N s + N d (number of ports) trajectories. The empirical observability Gramian requires N p + N q (number of states) trajectories. The empirical cross Gramian requires N s + N d + N p + N q trajectories, and the linear empirical cross Gramian requires 2(N s + N d ) trajectories. For the tested reductors this amounts to N s + N d trajectories for the POD, goal-oriented POD, and DMD-Galerkin method, while the port-Hamiltonian variants of the dominant subspaces, balanced POD, balanced truncation and balanced gains These predicted complexities are reflected in the offline runtimes, when computed sequentially. As the computation of trajectories is embarrassingly parallel, all trajectories are however computable simultaneously. Nonetheless, the complexities of the reachability-Gramian-only and port-Hamiltonian reductors are independent from the discretization, and thus most relevant for large-scale gas networks. Table 7 MORSCOREs μ(250, mach (16) ) in the L 2 ⊗ L 2 error norm for the "GasLib-134v2" benchmark network from Sect. 6

Outlook
The next stage in the development of morgen involves testing larger real-life networks, such as the deliverable of the SciGRID_gas 15 (Open Source Model of the European Gas Network) project. Yet, various further venues of linked modeling and model reduction questions are still not covered by morgen. In terms of model reduction, an interesting issue are intraday switchable valves, which change the topology of the gas network graph and likely require to extend the utilized model reduction methods towards these switched systems.
Another interesting question in need of further investigation is the minimal time horizon of the training phase. A lower bound is the time step times the longest path from a supply to a demand node, but this is likely not sufficient.
Besides an additional hyper-reduction module (Sect. 3.3) post-processing the reduced order models, a decoupler module pre-processing (DAE) models as described in [11,12] is projected.
Also as detailed in Sect. 3.4, the pipe roughness is a relevant parameter for (transient) simulations [130], yet the entailing high-dimensional parameter space, due to the locally differing roughness and attrition rates, would have to be treated, too. This in turn would raise the question for combined state and parameter reduction [65], and is postponed to future investigations.
Finally, using a tunable efficiency factor [105,107,111] that scales the model's friction term, can be used to tune the models to match real data.

Conclusions
In more than half a century of computational transient gas network simulation research and industrial use, morgen seems to be the first open-source platform covering modeling, simulation, and model order reduction of gas (and energy) networks. The target applications for morgen are finding the best model reduction method or best reduced order model for a network by heuristic comparison, as well as comparing model-solver-reductor simulation ensembles.
From a mathematical point of view, a next generation gas network simulation stack should consist of a (endpoint) port-Hamiltonian model, a (first-order) IMEX solver, and a block-diagonal Galerkin projection reductor, which is confirmed by the numerical results.
This results in the following heuristically determined but theoretically explainable recommended combination: The endpoint model together with the first-order IMEX solver, and a Galerkin reductor, specifically a structured dominant subspaces or structured DMD-Galerkin, exhibiting the highest MORscores in the numerical experiments. The performance of structured balanced truncation and the related structured balanced gains may be improved in terms of stability(-preservation) either by a variant of the technique [88], a stabilizing inner product [124], an (energy-)stable inner product [78], or an optimizationbased post-processing as in [21].
Lastly, we invite researchers, engineers and users to provide their reductors, solvers, networks and scenarios for expansion and testing with morgen for a broader view of this comparison. residence at the Institute for Computational and Experimental Research in Mathematics in Providence, RI, during the "Model and dimension reduction in uncertain and dynamic systems" program. Open Access funding enabled and organized by Projekt DEAL.

Availability of data and materials
The Matlab language source code of the morgen platform 1.0 is licensed under BSD-2-CLAUSE LICENSE, can be obtained from: https://doi.org/10.5281/zenodo.5012357 and is authored by: C. HIMPE and S. GRUNDEL.