{'fraction_non_alphanumeric': 0.05007185114681895, 'fraction_numerical': 0.01455685571916466, 'mean_word_length': 4.206373882627283, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On snapshot-based model reduction under compatibility conditions for a nonlinear flow problem on networks

This paper is on the construction of structure-preserving, online-efficient reduced models for the barotropic Euler equations with a friction term on networks. The nonlinear flow problem finds broad application in the context of gas distribution networks. We propose a snapshot-based reduction approach that consists of a mixed variational Galerkin approximation combined with quadrature-type complexity reduction. Its main feature is that certain compatibility conditions are assured during the training phase, which make our approach structure-preserving. The resulting reduced models are locally mass conservative and inherit an energy-bound and port-Hamiltonian structure. We also derive a well-posedness result for them. In the training phase, the compatibility conditions pose challenges, we face constrained data approximation problems as opposed to the unconstrained training problems in the conventional reduction methods. The training of our model order reduction consists of a principal component analysis under a compatibility constraint and, notably, yields reduced models that fulfill an optimality condition for the snapshot data. The training of our quadrature-type complexity reduction involves a semi-definite program with combinatorial aspects, which we approach by a greedy procedure.

reduction methods. As for the complexity reduction, we refer to the empirical interpolation method and its variants [6,18] as well as to the quadrature-type approaches [4,27].
While reduction methods have been successfully used in numerous applications, it is known that the standard approaches can suffer from poor results and stability issues for problems of complex structure. One major cause for the issues is that the reduced models may disregard fundamental structural properties of the model they approximate. These properties might include, e.g., conservation laws, dissipative relations or symplecticity. As a remedy, structure-preserving approximation methods have been developed, see, e.g., [1,15,19] for selected overviews. Problem-adapted space discretization approaches range from finite volume methods [44] to mixed finite element methods [10] or mimetic finite differences [28]. As for the time discretization, e.g., the symplectic and geometric integrators [33,41] are famous structure-preserving methods. The development of structure-preserving reduction methods is a more recent topic. Frequency-based model order reduction methods for linear port-Hamiltonian models are discussed in [32,53]. The works [13,14] treat nonlinear Lagrangian dynamics by snapshot-based approaches, particularly also using empirical interpolation-type complexity reduction. The papers [1,2,48] investigate model order reduction methods that preserve canonical symplectic structure under certain compatibility conditions and thus yield reduced systems in Hamiltonian form. Additionally, complexity reduction by the discrete empirical interpolation method is considered, but this step is not strictly structure-preserving, i.e., it does not guarantee a Hamiltonian structure. A modification of the discrete empirical interpolation that allows for a Hamiltonian representation can be found in [16], but it inherits higher local errors due to an enforced symmetrization step it includes. In other contexts [27,35], quadrature-type complexity reduction methods are used, which show to be structure-preserving in a more natural way.
In this paper, we construct structure-preserving and online-efficient reduced models for a nonlinear flow problem on networks. The problem is governed by the barotropic Euler equations with a friction term and inherits a port-Hamiltonian structure. It finds application in the context of gas transport networks. The efficient simulation, control and optimization of gas networks is an important topic as the wide range of publications in this direction shows. In [11] a model hierarchy of gradually simplified models is derived, whereby the isothermal (i.e., barotropic) Euler equations are the starting point. Strongly related are the contributions [20,21], which apply adaptive model-switching in a gas simulation tool. Other approaches to speed up calculations include problem-adapted preconditioners [49] or concepts from discrete optimization [40]. Conventional snapshot-based model reduction methods have also been applied for gas networks [31,38]. Our model reduction approach focuses on structure-preserving properties, i.e., local mass conservation, an energy bound and a port-Hamiltonian structure underlying the model problem. We establish these properties using compatibility conditions related to the mixed variational formulation [46], see also [22,24]. While beneficial for the robustness and performance of the reduced models, the compatibility conditions pose a challenge in the training phase, as one faces constraints in the training problems related to the snapshot data. We aim for reduced models that fulfill an optimality condition under these compatibility constraints. The conventional model reduction approaches, cf. [6,8,42], are designed towards an optimality condition related to the snapshot data, but compatibility conditions are not taken into account, i.e., structural properties are generally not guaranteed. For some applications, more problem adapted approaches have been derived. In [3,37,51] inf-sup stable reduced models for the Stokes and Navier-Stokes equations are derived. This is done using a splitted procedure that treats the training towards the snapshot data and the inf-sup compatibility separately from each other. Consequently, no optimality condition related to the snapshot data can be guaranteed in this approach. Structure-preserving model order reduction methods that aim for optimal reduced models have been considered for canonical Hamiltonian systems, see [1,48]. Conceptually, our approach shares most similarities with the latter approaches, which are also known as symplectic model reduction. However, the model problem we treat is quite different. Particularly, we cannot reformulate our problem as a canonical Hamiltonian system, we have to deal with network aspects, and our approximations rely on a mixed variational formulation. Regarding the model order reduction, we show that our training problem with its compatibility constraints can be attributed to an unconstrained principal component analysis, given appropriate norms are chosen. Based on this, we derive an efficient algorithmic implementation. Similarly, a training problem involving the snapshot data and constraints reflecting another compatibility condition is considered for the complexity reduction. It is a semi-definite program with combinatorial aspects, which we approach algorithmically by a greedy procedure. The performance of our proposed model reduction method is investigated numerically at the example of gas networks and compared to conventional, non-structure-preserving reduction methods. Moreover, we show a well-posedness result for our reduced models.
The structure of the paper is as follows. In Section 2 we state the model problem and our notation concerning the network aspects. The approximation ansatz from [46] is shortly reviewed in Section 3. It is employed to set up the underlying full order model (space discretization) and also serves as a basis for the reduction steps. Our structure-preserving model order-and the complexity-reduction methods are derived and analyzed in Section 4 and Section 5, respectively. The numerical studies are presented in Section 6, and a well-posedness result for our reduced models is derived in Section 7.

Model problem
This section introduces our nonlinear model problem, which describes, e.g., gas transportation networks [11,36]. Moreover, the variational principle our approximations are based on is stated, and the notation concerning the network aspects is settled, cf. [23,45].
Let a (gas pipe) network be described by a directed graph (N , E) with sets of nodes N = {ν 1 , . . . , ν } and edges E = {e 1 , . . . , e k } ⊂ N × N . Each edge e ∈ E is attached with a length l e > 0 and a weight A e > 0. The set of all edges adjacent to the node ν is denoted by E(ν) = {e ∈ E : e = (ν,ν), or e = (ν, ν)}, and a weighted incidence mapping is defined by A e for e = (ν,ν) for someν ∈ N −A e for e = (ν, ν) for someν ∈ N .
The nodes are grouped into interior nodes N 0 ⊂ N and boundary nodes N ∂ = N \N 0 . Function spaces on the network are constructed by compositions of standard Sobolev spaces for every edge. The spatial domain of the union of edges is Ω = {x : x ∈ e, for e ∈ E}. Every edge e can be identified with an interval (0, l e ) which is tacitly employed in the upcoming integral expressions. The space of square-integrable functions on E is given by L 2 (E) = b : Ω → R with b |e ∈ L 2 (e) for all e ∈ E , where the subscript . |e indicates the restriction of a function to the edge e. The respective scalar product and norm read b,b = e∈E A e e bb dx and ||b|| = b, b for b,b ∈ L 2 (E). The weak (broken) derivative operator for functions on the network is defined by The space of functions with square-integrable weak broken derivative is given as = 0, for ν ∈ N 0 } incorporates certain coupling conditions at inner nodes. The boundary nodes N ∂ = {ν 1 , . . . , ν p } are assumed to be connected to exactly one edge each. Thus, a boundary operator . . , p and b ∈ H 1 pw (E). We consider the barotropic Euler equations with friction governing density and velocity ρ, v : [0, T ] × Ω → R on the edges of the network, Here, P is the derivative of the pressure potential P : R + → R, which we assume to be strictly convex and twice continuously differentiable for all considered ρ. Moreover, the friction termr is assumed to be pointwise positive. The solution components are interconnected by the coupling conditions at ν ∈ N 0 The model is closed by appropriate boundary and initial conditions, e.g., and . Instead of the pressure potential, the flow can also be described by means of the pressure p given by p(ρ) = ρP (ρ) − P (ρ), cf., [5,29]. The physical energy of the system is given by the HamiltonianH Conservation of mass as well as of the Hamiltonian at inner nodes are ensured by the coupling conditions. Furthermore, for sufficiently smooth solutions, it can be shown that the Hamiltonian dissipates over time, up to the exchange with the boundary, i.e., The energy dissipation (in-)equality can be derived from the variational principle stated in (3). For a strong solution (ρ, v) ∈ C 1 ([0, T ]; C 1 pw (E) × C 1 pw (E)) of (1) and the mass flux m = ρv, the variational principle holds for all q ∈ L 2 (E), w ∈ H 1 div (E), whereby u = [u ν1 , . . . , u νp ] T and r is given by r(ρ, m) =r(ρ, m/ρ). The principle follows from (1) by standard arguments and the coupling conditions inherited in H 1 div (E). Remark 2.1. The dissipation (in-)equality (2) is a consequence of the more profound structural property that the barotropic Euler equations can be stated in port-Hamiltonian form. In particular, (1a) can be formally written as with ∇ ρH and ∇ vH denoting the functional derivatives ofH with respect to (w.r.t.) the variable ρ and v, respectively. The approximation procedure from [46], which we employ in this work, is designed to preserve the port-Hamiltonian structure on the discrete level. It crucially relies on the formulation of the variational principle (3) in terms of the mass flux m while still using the Hamiltonian formulation with its symmetries, as it is stated here.

Structure-preserving approximation framework
The structure-preserving approximation concept from [45,46] forms the theoretical basis of this work. We briefly present the concept and then discuss our specific choices for constructing the full order model and setting up the data-based reduction.

Generic approximation procedure
The structure-preserving approximation procedure developed in [45,46] is based on the variational formulation (3). It is applicable to Galerkin projections with classical spatial finite element discretization and model order reduction as well as to complexity reduction of quadrature-type for the nonlinear terms. On all approximation levels local mass conservation and an energy bound for the approximations as well as the port-Hamiltonian structure indicated in Remark 2.1 are kept under mild assumptions.
The space K consists of the edge-wise constant functions fulfilling a certain coupling condition. Thus it only depends on the underlying topology and dim(K) = |E| − |N 0 |. Assumption 3.2 (Compatibility of scalar product). Let the bilinear form ·, · * : L 2 (E) × L 2 (E) → R be such that the following holds: for all q ∈ Q, w ∈ W. Further, H * (ρ, m) = P (ρ)+m 2 /(2ρ), 1 * is referred to as the Hamiltonian of the system. Theorem 3.4 (Energy bound, [46]). The energy dissipation (in-)equality is fulfilled by any solution (ρ, m) of System 3.3, Note that the inequality in the energy bound is implied by the non-negativity of the map f → f, 1 * . Our flow problem possesses the form of System 3.3 on each approximation level. At first, a Galerkin projection is constructed by mixed finite elements yielding the full order model (FOM). Subsequently, a second Galerkin projection is applied to obtain a reduced order model (ROM). Although the ROM is of much lower dimension than the FOM, it is unfortunately not more efficient, since the evaluation of the nonlinearities is not independent of the dimension of the FOM. Thus, a final complexity reduction step is performed that replaces the L 2 -scalar product by a cheaper-to-evaluate bilinear form for all nonlinear expressions yielding the complexity-reduced model (CROM). In this sense, ·, · * in System 3.3 is chosen as a complexity-reduced bilinear form (specified and denoted as ·, · c below) for the CROM, otherwise for FOM and ROM it is considered as the L 2 -scalar product ·, · . In any case, ·, · * is a scalar product and induces a norm on the spaces Q and W, respectively, by Assumption 3.2. The energy dissipation (in-)equality (Theorem 3.4) ensures the stability of the model on each approximation level. Apart from that, System 3.3 inherits further structural properties. It preserves the port-Hamiltonian structure mentioned in Remark 2.1, and its solutions are locally mass conservative, i.e., ∂ t ρ = −∂ x m holds in a pointwise sense. We refer to [45,46] for details and proofs.

Specific full order model and data-based reduction setting
Concerning the FOM we follow [24,25] and use a mixed finite element space that is compatible in the sense of Assumption 3.1. A partitioning of the spatial domain Ω into the cells K 1 , . . . , K J is defined such that j K j = Ω, Kj ∩K l 1dx = 0, for j = l, K j ⊂ e for a e ∈ E for j = 1, . . . , J. (4a) The finite element ansatz space V f = Q f × W f , which we also refer to as FOM space, is given by The functions in Q f are piecewise constant and discontinuous, whereas W f consists of piecewise linear, edgewise continuous functions and inherits the coupling conditions that relate to mass conservation at inner nodes. The reduced models are trained towards data, the so-called snapshots, given as The snapshots are typically selected (time) points of one or several solution trajectories of the FOM for appropriate training scenarios. In the model order reduction we seek for a low-dimensional subspace dim(V f ), that fulfills Assumption 3.1 and in which the snapshot data is approximated with high fidelity. The subsequent complexity reduction aims for constructing a bilinear form ·, · c that fulfills Assumption 3.2, whose evaluation-cost is independent of dim(V f ), and for which f (ā ), w r c ≈ f (ā ), w r holds for = 1, . . . , L and w r ∈ W r for the appearing nonlinear terms f . We use a quadrature-type approximation.
Our approach differs from most conventional data-based model reduction methods, e.g, [8,18,43], in that compatibility conditions are imposed to ensure the preservation of structural properties for the reduced systems, such as an energy bound (Theorem 3.4) or well-posedness (Section 7). In the training phase, the conditions pose an additional challenge, as they lead to constraints.

Realization of model order reduction
In this section we formulate the data-based model reduction problem with compatibility conditions. It can be considered as a principal component analysis (PCA) with constraints. As will be shown, this constrained problem can be attributed to an underlying standard PCA given suitable norms are chosen. Based on that, we derive an efficient algorithmic implementation.

Constrained principal component analysis
Our model order reduction follows the basic notion of proper orthogonal decomposition. The ROM space V r is obtained by a PCA of the snapshot data, but in contrast to the standard approach, constraints related to the compatibility conditions are posed. To be able to treat the constraints efficiently, the choice of a proper scalar product ·, · for W f in the PCA will be crucial. We denote the respective induced norm by || · || , i.e., ||w|| = w, w . Moreover, the orthogonal projection onto a subspace U w.r.t. the scalar product ·, · and ·, · is denoted by Π U and Π U , respectively. The general approximation problem we consider for the determination of the ROM space V r reads as follows.
Due to the compatibility conditions (Assumption 3.1), Problem 4.1 is a constrained optimization problem which might be quite involved to solve. We aim for the derivation of an equivalent unconstrained standard PCA. To remove the constraints, the underlying pure data approximation problem has to be separated from the question of compatibility.
, the assertion is straightforward to show. We refer to [23], where the norm || · || is considered.
In view of Assumption 3.1, we introduce the orthogonal decomposition of W f w.r.t. ·, · into K and its orthogonal complement K ⊥ , Furthermore, we define ∂ + x : Q f → W f as the right-inverse of the spatial derivative ∂ x : W f → Q f , which additionally fulfills Π K ∂ + x ≡ 0. Note that K is the kernel of ∂ x , and that ∂ + x is well-defined and relates to a weighted Moore-Penrose inverse with weights given by the scalar product ·, · , cf. [7]. First, we show that the constraints can be formally removed for any choice of ·, · .
be a compatible reduced space in the sense of Assumption 3.1. Analogously as for the FOM, an orthogonal decomposition of W w.r.t. ·, · can be performed, which yields The orthogonal decomposition of the FOM implies ||w|| 2 = ||Π K ⊥ w|| 2 + ||Π K w|| 2 for any w ∈ W f . Using this equality for the approximation error together with the compatibility All in all, the compatibility conditions of Assumption 3.1 can be removed, and the unconstrained optimization problem is equivalent to Problem 4.1 in the sense of the theorem.
Among others, the theorem shows that the snapshots m should not be used directly in the PCA but should be altered first to Π K ⊥ m . Moreover note that the sub-problem determining Q r is unconstrained. However, its cost functional inherits terms with different norms depending on the choice of scalar product ·, · . For the specific choice (6), we can attribute the following standard PCA to the training problem.  (5) and Theorem 4.5. Let ·, · be the scalar product defined in (6).
Proof. As auxiliary result we show that for any compatible pair of spaces Q, W it holds By the definition of orthogonal projections, The last equivalence makes use of K being the kernel of ∂ x . The full equivalence together with the fact ∂ x w 2 ∈ Q yields that ∂ x w 2 is the orthogonal projection of ∂ x w onto Q, i.e., the result (7).
To show the equivalence of the minimization problems, we make use of the representation given in Theorem 4.3. Its cost function includes terms of the form For our specific choice of inner product, it holds whereby (7) has been used in the last step. From that the equivalence of the cost function of Problem 4.4 and the one stated in Theorem 4.3 can be concluded in a straightforward manner. Using Theorem 4.3 finishes the proof.

Structured representations
For the algorithmic realization we introduce coordinate representations of the full and reduced order models.

Full order model
The mass matrices for Q f and W f as well as the coordinate representation of ∂ x : W f → Q f are defined by The boundary operator T : W f → R p has the coordinate representation B = T w 1 , . . . , T w N2 T . Moreover, given (ρ, m) = Ψ(a) ∈ V f , we define vectors related to the nonlinear expressions by Then the FOM (System 3.3) can be equivalently described by the following algebraic representation.
System 4.6 (Algebraic representation of FOM). Given appropriate a 0 ∈ R N and u : Remark 4.7. Related to System 4.6, the Hamiltonian H can be defined as H(a) = H(Ψ(a)). Using this Hamiltonian, the system can be transformed into standard port-Hamiltonian form, see [46].
which can be concluded from Assumption 3.1 (∂ x W f = Q f ).

Reduced order models
The projection of the FOM onto the low-dimensional ROM space V r ⊂ V f can be equivalently described in the algebraic setting as projecting System 4.6 with a reduction matrix V ∈ R N,n and n = dim(V r ) N . As we consider projections orthogonal w.r.t. the L 2 -inner product, weighted orthogonal projections occur naturally in the coordinate representation. Accordingly the following definitions are helpful, cf., [7]. Given a symmetric positive definite matrix M ∈ R N,N , we introduce the scalar product ·, · M and the respective norm · M , i.e., x, y M = x T My and x M = x, x M for any x, y ∈ R N . For any matrix V ∈ R N,n of full column rank, we refer to as the weighted orthogonal projection w.
The compatibility condition of Assumption 3.1 on V r can be recast as a compatibility condition on the reduction matrix V according to Lemma 4.8, as can be derived with the help of (10), cf., [25,45]. The ROM and a reduced coordinate representation a r , Va r ≈ a, are then characterized by a projected version of System 4.6. Its initial conditions a r (0) ∈ R n1+n2 are chosen as the L 2 -projections of a 0 = [a T 0,1 , a T 0,2 ] T ∈ R N1+N2 , accordingly. System 4.9 (Algebraic representation of ROM). Given System 4.6 and a reduction matrix V as in Lemma 4.8, where the reduced state matrices are

Computation of structure-preserving reduction basis
To reformulate a PCA problem in a weighted norm in terms of the Euclidean norm, we make use of the following result, cf. [26]. Let M ∈ R N,N be a symmetric positive definite matrix and M = L T L. Then it holds The relation follows from . Then Problem 4.4 can be equivalently described by where the vectorsṽ i denote the first left-singular vectors of the matrix K.
Proof. The stated minimization problem relates to a coordinate representation of Problem 4.4 in the basis {q 1 , . . . , q N1 } of Q f . This can be seen from the equalities according to (12). Hence, V * 1 solves the problem of the theorem, if and only ifṼ * The latter is a standard PCA problem and can be solved by the method of snapshots, cf. [42,43]. It has a solutioñ V * 1 = [ṽ 1 , . . . ,ṽ n1 ] withṽ i as in the theorem. Consequently, V * 1 fulfills im(V * 1 ) = im(L −1 [ṽ 1 , . . . ,ṽ n1 ]). In the following algorithm we summarize our procedure for constructing a reduction matrix V, which yields an optimal reduction space in the sense of Problem 4.1. As the characterizations of V r = Q r × W r in Lemma 4.8 and Theorem 4.3 suggest, the procedure separates into two parts. First, the space Q r = im(V 1 ) is determined according to Theorem 4.10. Then, W r = im(V 2 ) is derived using the compatibility conditions. (3) Determine matrices V i ∈ R Ni,ni for i = 1, 2 (whose column vectors are orthonormal in the scalar products induced by Q and W, respectively) with Step (3) of Algorithm 4.11, weighted scalar products are taken into account for better numerical stability. Particularly, J † is the right-inverse of J weighted w.r.t. the scalar product induced by W. Let us stress that the equality im(V 2 ) = ker(J) ⊕ im(J † (QV 1 )) as well as im(V 1 ) and im(V 2 ) are independent of the choice of scalar product in the absence of rounding errors.
Remark 4.12. The computational complexity of Algorithm 4.11 is comparable to the standard POD approach when implemented appropriately. Since Q is diagonal and W is tridiagonal, their Cholesky factorizations are of the same form.
Step (2) scales only linearly in the dimension of the FOM when the method of snapshots [42,43] is employed. In Step (3), the kernel ker(J) can be efficiently determined by a sparse LU-factorization [30]. Moreover, the right-inverse J † does not have to be computed explicitly, but a few sparse linear equations can be solved instead, as can be seen using the Cholesky factorization of W.

Quadrature-type complexity reduction
Although the ROMs are typically of much lower dimension than the FOM, i.e., dim(V r ) N = dim(V f ), they are in general not more efficient due to the lifting bottleneck related to the nonlinearities. The evaluation of the nonlinear integrals in System 3.3 scales with the number of finite elements K j , j = 1, . . . , J, and hence is not independent of the FOM dimension, since N = 2J + dim(K), see, e.g., the friction term We propose a complexity reduction by a quadrature-type approximation that restricts to integral evaluations at only few finite elements. It is trained towards the given snapshot data and, in contrast to conventional methods, regards compatibility conditions (Assumption 3.2). The resulting training problem is a semi-definite program with combinatorial aspects. We approach it algorithmically by a greedy procedure, as is generally done for complexity reduction.
Remark 5.1. The CROM distinguishes from the ROM only in the nonlinear terms. Its algebraic representation is similar to System 4.9, but f α r , f β r , f γ r are substituted by complexity-reduced analogs.

Training goal
The complexity reduction is done subsequently to the model order reduction and its training employs the ROM space V r = Q r ×W r explicitly. Similarly to [27,35], we consider a complexity reduction of quadrature-type. We aim for a complexity-reduced bilinear form that approximates ·, · , satisfies Assumption 3.2 and realizes ·, · * in System 3.3. In the training phase, we search for the index-set I ⊂ {1, . . . , J} and weights ξ i ≥ 0 for i ∈ I. In particular, we aim for |I| = n c with n c J small to ensure a reduced computational complexity for the CROM. To formulate the training goal, we assume an appropriate collection of snapshots of all nonlinear integrands in the ROM to be given, i.e., The complexity reduction is desired to be of high fidelity for the snapshot data, i.e., f , 1 c ≈ f , 1 for = 1, . . . ,L. Expressing that in a least-squares sense yields the following optimization problem for the training of the quadrature-rule.
Problem 5.2. Let nonlinear snapshot functionsS f be given as in (14). Let n c ∈ {1, . . . , J} andC > 1. Find an index-set I * ⊂ {1, . . . , J} and weights ξ * i ∈ R for i ∈ I * as solution to the problem min I⊂{1,...,J}, |I|=nc ξi≥0, i∈IL =1 i∈I The constraints reflect the compatibility condition of Assumption 3.2 that ensures structure preservation. In particular, the non-negativity of ξ i directly implies the non-negativity of the function f → f, 1 c .
As for the snapshotsS f , we consider the following mappings forā = (ρ, m), related to the three nonlinear integral expressions occurring in the ROM (System 3.3), with w j r for j = 1, . . . , n 2 denoting the basis for W r . For the first two mappings, nonlinear snapshots are generated using the data set of the state snapshotsS = {ā 1 , ...,ā L },ā ∈ V f , i.e., the snapshots the ROM is trained for. The third mapping appears in a time-derivative of the ROM. To treat this accordingly, we additionally assume that the state snapshotsS consist of one or several solution trajectories in time. Particularly, givenā l andā are snapshots of the same trajectory at different time points t k−1 and t k , we collect the terms (v(ā ) − v(ā l ))/(t k − t k−1 )w j r for j = 1, . . . , n 2 , which relates to using a finite difference approximation on ∂ t v(ā) at time t k . In case thatS represents one trajectory, we can approximate the time-derivative at L − 1 instances, and the collection of all nonlinear snapshotsS f hasL = (3L − 1)n 2 entries in total. Remark 5.3. By construction, the number of integrand-snapshots inS f is significantly larger than the number of state snapshotsS in (5). When the problem at hand is of very large scale, a pre-processing ofS f or other strategies to speed up the training phase might become necessary. We do not discuss this issue here but refer to [4,35] for some attempts in this direction.

Greedy implementation of training
For algorithmic reasons, we rewrite Problem 5.2 in algebraic form. Let ξ = [ξ 1 , . . . , ξ J ] T ∈ R J denote the extended vector of weights, which is composed of the weights ξ i for i ∈ I used in the quadrature and the entries ξ j = 0 for j / ∈ I. With f ∈ R J defined by [f ] i = Ki f dx for the nonlinear integrands f ∈ L 1 (E), = 1, ...,L in (14), we can rewrite the cost functional in terms of ξ and The constraint in Problem 5.2 enforces a norm equivalence between the L 2 -norm and its complexity-reduced counterpart || · || c . And as the ROM bases {q 1 r , . . . , q n1 r } and {w 1 r , . . . , w n2 r } are orthonormal in the L 2 -norm according to Algorithm 4.11, the related mass matrix becomes the unit matrix. Its complexity-reduced counterpart reads The constraint (17c) makes this a semi-definite program [54]. And inherently in the complexity reduction, there is a combinatorial aspect, in our case the selection of the index set I of active finite elements. Particularly the latter makes this problem too hard to be solved to global optimality. Thus, we rely on a greedy approximation, as most complexity-reduction methods do, cf., [2,6,18], see also [4,27] for quadrature-type approaches similar to ours. The basic idea is to alternate between enlarging the index set I by a greedy search and constructing an optimal weight for that fixed index set. Algorithm 5.4 summarizes the procedure. While the non-negativity constraint for the weights ξ i is included explicitly, the eigenvalue constraint (17c) is only included as a safeguard at the end of the algorithm. As will be shown in the numerical results in Section 6, the latter condition is obtained by our greedy search for reasonable choices of n c without further doing for our problem at hand.
Set n c ← n c + 1 and go to (3a) endif end Let us note that only n c columns or less of the data matrix A are used in most steps of the algorithm. The exception is Step (3c), where the gradient is evaluated once. It can be avoided to store the full matrix, if storage requirements are critical, cf., [9,47] and Remark 5.3. Most importantly, the non-negative least-squares problem of Step (3d) needs only to be solved over very low-dimensional linear spaces. We do this using the built-in Matlab routine lsqnonneg in our implementation, which realizes an active-set method.

Numerical results
For the numerical studies of this section we consider our model problem in parameter regimes relevant for gas distribution networks. The model equations are given by the isothermal Euler equations in a friction-dominated regime. The network parameters are extracted from the gaslib [52] with minor adaptions. Networks of pipes are considered, where each edge e ∈ E represents a pipe of a cross-sectional area A e = π/4 (D e ) 2 . The General boundary conditions are treated according to [46]. Specifically, for each mass flow boundary condition a Lagrange multiplier is added and also kept in the reduced models. The FOM is realized according to Section 3.2 using a uniform mesh on each pipe with the maximal spatial step size ∆ x ≤ 200 [m]. Time discretization is carried out by an implicit Euler method in the primitive variables ρ and v, cf. [46], with a constant time step ∆ t = 1 [s], i.e., 3600 steps per hour of simulation time. The resulting nonlinear systems in each time step are solved by the Newton's method. If this fixed-point iteration diverges for any time step, we consider this as a simulation-breakdown. Given the FOM solutionā, the error of an approximationã (obtained by a reduced model or orthogonal projection) is measured by the L 2 -norm in space and the supremum-norm in time, yielding the relative error E T = (max t∈[0,T ] ||ā(t) −ã(t)||)/max t∈[0,T ] ||ā(t)||. All numerical results have been generated using MATLAB Version 9.1.0 (R2016b) on an Intel Core i5-7500 CPU with 16.0GB RAM.
In Section 6.1 some qualitative properties of our complexity reduction approach are showcased at a smaller academic network scenario (diamond network). More quantitative studies with comparisons to other wellestablished, but non-structure-preserving, model reduction methods are presented in Section 6.2 using a network of realistic size.

Qualitative study for complexity reduction
The topology and parameters for the academic network example (diamond network) are given in Fig. 1. For the friction factor, we consider two different choices λ = 0.01 (high friction) and λ = 0.002 (lower friction). At the two boundary nodes ν 1 , ν 2 , we prescribe the boundary conditions according to with reference time t * = 1[h] and the time-varying input profile The initial values are chosen as the stationary solution related to the boundary conditions at t = 0. The FOM yields a system of dimension N = 1922, and the reduction methods (model order-and complexity-reduction) are trained towards 1000 equally distributed time snapshots of the solution trajectory of the FOM. The dimension of the ROM space V r is fixed to n = 14. As regards the complexity reduction, we vary the parameter n c , as the qualitative behavior of this step is is the focus in the upcoming. For the two different cases of friction factor, the time evolution of the solution is visualized in Fig. 2. Hereby, the end of the edge e 2 (at node ν 4 ) is used as fixed spatial position. As to be expected, stronger damping effects are observed for larger λ, especially for the density. Generally speaking, models with higher damping effects are mostly better suited for model reduction, but the effect seems rather minor in the parameter range we encounter in the context of gas transportation. For both choices of λ, the CROMs with n c = 20 produce very accurate time responses, see Fig. 2. But also when the complexity reduction is not well-resolved, as is the case for n c = 14, stable simulations are produced, which is in accordance to our energy bound (Theorem 3.4). The relative errors over the full simulation with λ = 0.002 are shown Fig. 3-right. We observe that the error the complexity reduction step adds becomes negligible for n c ≥ 18. Moreover, the compatibility condition on the spectrum of the complexity reduced mass matrix, σ(M c ) ⊂ [C −2 ,C 2 ], is fulfilled for n c ≥ 12 and moderately largeC, as is indicated by the condition number of M c , see Table 1.
An important qualitative property concerning the training of our complexity reduction method, which we want to showcase, is its adaption to the ROM space. To demonstrate that, we consider reduced models with an artificially enlarged ROM space given asṼ r = V r +Ũ, whereŨ consists of functions with very local support. That is, two of the basis functions of the FOM for the density, which lie symmetrically around the midpoint  of edge e 3 , are added to Q r , and the space W r is supplemented according to Assumption 3.1-(A1). ThusṼ r is compatible, and it holds 18 = dim(Ṽ r ) > dim(V r ) = 14, but almost no improvement in terms of fidelity is observed for the scenario. In fact, the extensionŨ causes the complexity reduction to produce less accurate models with larger cond(M c ), at least for small n c , see Fig. 3-right and Table 1. The reason becomes evident when looking at the quadrature points obtained for the CROM without and withŨ. In Fig. 3-left the locations and values of the quadrature weights are visualized along the path of the edges e 1 to e 4 (marked in magenta in Fig. 1). When the spaceŨ is added, we find three weights located on its support, see the encircled markers  in the figure. Let us also mention that the weights do not cluster, the pairs that are observed around the junctions (the cyan vertical lines) lie on different pipes. Summarizing, the enrichment of the reduced space bỹ U hampers the complexity reduction, but our method still produces stable simulations. This illustrates that our training procedure promotes Assumption 3.2 and can cope with an enrichment of the reduction basis with elements chosen independently of the snapshot data. In more relevant scenarios, the space enrichment might be necessary due to compatibility conditions such as Assumption 3.1, and the effects on the complexity reduction can be expected to be less exaggerated. Let us stress that other complexity reduction methods, such as the empirical interpolation methods [6,18], are trained independently of the choice of the ROM space and thus cannot adapt to it. Possibly related to that, we observe that the discrete interpolation method (DEIM) inherits stability issues when combined with our model order reduction approach, see the end of Section 6.2.

Comparison to non-structure-preserving reduction methods
The benefits of our structure-preserving approach are demonstrated for a benchmark with a realistic gas network, illustrated in Fig. 4. The network is a slight modification of [52,, cf. [46]. Let us note that a few compressors and valves would be needed in an actual application. We do not consider these active elements for ease of presentation, as they are prescribed by low-order algebraic models, which would not be regarded during the reduction process anyway. The network consists of 38 pipes with diameters D e between 0.4 and 1 [m] and lengths l e between 5 to 74 [km]. The total pipe length is 1008 [km]. Over the whole network, we set the friction factor to λ = 0.008. At the six boundary nodes ν i , i = 1, ..., 6, we prescribe the following boundary conditions for t ∈ [0, 5t * ] and reference time t * = 1[h], ρ(t, ν 1 ) = (65 + u(t/t * )) ρ , ρ(t, ν 2 ) = (50 + u(t/t * )) ρ , ρ(t, ν 4 ) = (60 − u(t/t * )) ρ , ρ(t, ν 5 ) = 60 ρ , ρ(t, ν 6 ) = 45 ρ , Am(t, ν 3 ) = −100 (Am) , with the input profile u ∈ {u A , u B } varied over two cases according to see Fig. 5. As initial condition the stationary solution belonging to the boundary conditions at t = 0 is taken. The FOM is a system of dimension N = 10156. The training of all reduction methods is realized with the trajectory of the FOM solution for Case A with 1000 equally distributed snapshots. Thus, we refer to Case A as the (perfectly) trained case and to Case B as the not trained one. In general, one can expect model reduction methods to perform worse in a not trained case, even more so if they do not have good stability qualities. In the upcoming, we first show temporal and spatial visualizations of the reference solution, which is given by the FOM. Then we compare our structure-preserving model order reduction approach against a conventional method, and subsequently do the same for the complexity reduction step. The time response of the FOM at the ends of pipe 2 and pipe 5 is shown in Fig. 6. Similarly as for the academic test scenario, we observe notable damping effects. The time response is smoothed out stronger at pipe-end 5 than at pipe-end 2, because the former has a larger distance to the driving boundary nodes. In Fig. 7, a spatial representation of the solution along the path from pipe 2 to pipe 8 (marked in magenta in Fig. 4) is shown for different time instances. The mass flows show to be discontinuous over junctions with more than two pipes meeting, as to be expected from the model. Apart from that, the spatial representations are rather smooth, which indicates that the solution manifold can be approximated well in low dimensions and model reduction methods could potentially be efficient.
We compare our proposed model order reduction against a conventional block-structured POD, which is very similar to the approaches [31,38] used in the gas network context. The reduction basis V BP for the conventional method is obtained by separately applying a principal component analysis to the density-and mass flux-snapshots and extracting reduction spaces of the same dimension for them (n 1 = n 2 ). This approach is preferable to applying POD onto the full state, but it is not structure-preserving in the sense of Assumption 3.1. In order to investigate the fidelity of the reduction spaces, independent of the stability properties of the reduced models, we additionally consider the L 2 -projections of the FOM solutions onto the ROM spaces. The respective projectors are denoted by Π Vr and Π V BP , and they yield the pure projection errors. The relative errors are shown in Fig. 8. There is one scenario, where the conventional POD shows slightly better results than our method. That is, better projection errors are observed for the perfectly trained case for the reduced dimension n small (Fig. 8-left, n ≤ 50). This is not surprising, as taking the L 2 -error favors V BP . Recall that our method is derived using another norm and regards compatibility conditions. Despite that, our reduction basis V r has slightly better projection errors in the trained case when the dimension is chosen to be larger, and shows significantly smaller projection errors for all parameter choices in the not trained case. The theoretical advantage our approach has over the other is that it fulfills an optimality condition in the full state. We suppose this is also the reason for the comparably low projection errors we observe. Practically more relevant is the comparison of the reduction errors. In this respect, our method outperforms the conventional reduction method much clearer, and we observe better errors for all dimensions n. The difference is most evident for the non-trained case, where the second smallest reduced order model with V r is of higher fidelity than any of the reduced models for V BP , see Fig. 8-right. Further, we want to highlight the small gap between reductionand projection-error for our proposed method compared to the conventional one, which strongly indicates the superior stability and robustness of our approach.  As regards the comparisons for the complexity reduction step, the starting point is a ROM obtained by our proposed model order reduction method. Our quadrature-type complexity reduction is compared to a nonstructure preserving alternative, which we representatively choose as the discrete empirical interpolation method    (DEIM) [17,18]. For convenience, we denote the DEIM-space dimension used in the training as n c . Note that the parameter in our CROM, which we also denote by n c , has a quite different interpretation, which is the number of nonlinear integrals that need to be evaluated. The DEIM is not directly related to the integral expressions, but rather to the algebraic representation of the system, and needs several nonlinear integral evaluation for each n c , which makes it less online-efficient in our setting. We apply the DEIM to each nonlinearity separately, i.e, given the algebraic representation of the ROM (System 4.9), that means the terms f α r , f γ r andf β r are separately complexity reduced. The latter termf β r is defined by f β r = J Tf β r , and it is introduced as the resulting DEIM approach regards some of the structural properties of our model problem, cf., Remark 6.1 and [45,46]. But let us emphasize that this DEIM approach still cannot guarantee port-Hamiltonian structure or a provable energy bound, in contrast to our quadrature-based approach. We test the complexity reduction methods for ROMs of dimensions n = 48 and n = 66, the resulting reduction errors are shown in Fig. 9. The first observation to be made is that DEIM yields in several parameter settings unstable results. Particularly in the untrained case, all DEIM models with n c ≤ 92 have either a simulation breakdown or very poor fidelity. Our CROM method shows to be much more robust and no simulation breakdowns occur. Let us also emphasize that we have a strong a priori indicator for stability for the CROM, which is the compatibility condition σ(M c ) ⊂ [C −2 ,C 2 ] for a reasonably small constantC. As Table 2 shows, the condition is fulfilled in our experiments. Another difference of the complexity reduction methods we want to stress is that our quadrature-based complexity reduction in CROM depends stronger on the underlying ROM dimension n. Particularly, when perfectly trained and n is chosen large, the error of DEIM becomes negligible for n c ≈ 92, whereas this holds only from n c ≈ 120 on for the CROM, see Fig. 9-left with n = 66. However, this does not mean that DEIM is more efficient in this or any other setting. As mentioned above, more nonlinear integrals for each n c have to be evaluated in DEIM than for CROM. This is directly reflected in the simulation times, which we can be found in Table 3 for the untrained case. The DEIM models have runtimes between 432 to 482 seconds, as compared to the CROMs with about 126 to 152 seconds. Thus, the CROM is more than three times faster than DEIM for models of similar fidelity. Compared to the 2253.6 seconds runtime of FOM, CROM shows a speedup of about 16.
To summarize, our quadrature-based reduction yields significantly more robust and efficient results as DEIM when combined with our model order reduction approach. Moreover, the compatibility condition of Assumption 3.2 is conveniently promoted by our greedy training procedure in all our tests. Remark 6.1. Instead of applying DEIM tof β as we did, one could apply it directly to f β . This alternative showed severe stability issues in our numerical tests. We assume this is related to a more profound loss of structural properties. Specifically, the anti-symmetry (symplectic structure) revolving around the terms −J and J T is destroyed by the more naive DEIM version, cf., [12,48] for related discussions.

Well-posedness of approximations
In this section, we derive a well-posedness result for our model order-and complexity-reduced approximation (CROM). The CROM is a realization of System 3.3, and its ansatz space V r is a subspace of a finite element space V f with partitioning K 1 , . . . , K J given as in (4). Further, the complexity-reduced bilinear form and Hamiltonian take, according to (13), the form with ξ i > 0 for i ∈ I, and H c (ρ, m) = h(ρ, m), 1 c .
Here, we assume ξ i > 0 instead of ξ i ≥ 0 for technical reasons, which can be done w.l.o.g. by simply restricting I accordingly. The well-posedness result is established under the situation described in System 3.3 (i.e., for boundary conditions as in (1c)) and the following additional assumptions on the model problem.
Assumptions (A1)-(A2) are used to show that H c is well-defined and bounds the norm of the solution. Moreover, it should be mentioned that P is also bounded from below by (A1). To establish a global existence result, we derive uniform boundedness of the solution and strict positivity of the density. The positivity of the density is needed to avoid zeros in the denominators of System 3.3. Our proof of that relies on (A2)-(A3). Further note that the restriction of the density to a single element K i , ρ(t) |Ki , is constant.
Remark 7.2. Similar results as the ones derived in this section can be established using other boundary conditions or friction/dissipation terms, but the derivations might get slightly more technical [45]. We also refer to [22,39] for related results for other space discretizations of the Euler equations. However, the latter references do not treat model order-and complexity-reduction. Lemma 7.3. Under the assumptions of this section, the solutions of System 3.3 fulfill ρ(t) |Ki > 0 for t > 0 and i ∈ I. Moreover, there exists a constant C, independent of the discretization parameters, such that where ∆ xi is the grid size of the finite element K i .
Proof. At several instances, we employ that ρ(t) |Ki and ∂ x m |Ki (t) are constant in space. As long as 1/ρ(t) |Ki is well-defined, it therefore holds where the subscript K i indicates a restriction of the spatial domain to one finite element. By the inverse estimate, there exists a constant C with ||∂ x m(t)|| Ki ≤ C/∆ xi ||m(t)|| Ki . Together with 1/ρ(t) and ∂ x m(t) both being constant on K i , this yields Setting together the two estimates, we thus get by the Gronwall lemma 1 ρ(t) |Ki ≤ exp t 0 c i (s)ds 1 ρ(0) |Ki .

It remains to bound
t 0 c i (s)ds for i ∈ I. As a preparation, we introduce the auxiliary functionĤ(ρ, m) = H c (ρ, m) + max 0 , −inf y∈(0,M ) P (y) , which is a constant shift of the Hamiltonian H c . By construction, it yields only non-negative values for all considered states and fulfills the same energy-dissipation equality as H c , cf., Theorem 3.4. Assuming ρ(t) |Ki > 0 and i ∈ I, we can follow wit R defined as in the lemma. The latter equality follows from integrating the energy dissipation equality for H (Theorem 3.4) in time. Moreover, R(t) ≥ 0 holds for t ≥ 0, as the last equality also implies that R can be bounded from below byĤ, which itself is non-negative. Now first applying the Jensen-inequality [50] and then inserting the former estimate yields Inserting the latter bound on t 0 c i (s)ds into the estimate obtained by the Gronwall lemma finishes the proof.
Remark 7.4. In the special case of ·, · c chosen as the L 2 -scalar product, we recover the case of pure Galerkin approximation, without complexity reduction. Lemma 7.3 then shows strict positivity of ρ(t) on all of the spatial domain. Such a result has been derived in [22] in a similar setting.
Note that positivity of ρ(t) |Ki can only be guaranteed for i ∈ I by Lemma 7.3. This turns out to be sufficient, as the ρ(t)-terms in the denominator are only evaluated for i ∈ I. Next, we derive a boundedness result for the solution.
Theorem 7.5. Under the assumptions of this section, there exist constants C 1 , C 2 , independent of the discretization parameters, such that for (ρ, m) ∈ V r with ρ |Ki > 0 for i ∈ I it holds Proof. Let k ∈ argmax i∈I ρ |Ki , and letĪ = {i ∈ I : ρ |Ki > 1]} w.l.o.g. be non-empty. Otherwise we trivially can bound the terms ||ρ|| c , ρ |K k by a constant. From Assumption 3.2 we get By Assumption 7.1-(A1) it follows that the function P can be bounded from below (not necessarily by zero), which implies that there exists a constantĉ such that ξ k P (ρ |K k ) ≤ i∈Ī ξ i P (ρ |Ki ) ≤ P (ρ), 1 c +ĉ.
Similarly, it follows where in the last step the same estimate on P (ρ |K k ) as before has been used. With the help of Young's inequality, it follows ||m|| ≤C √ ξ k ( P (ρ), 1 c +ĉ) Setting together the estimates for ||ρ|| and ||m|| shows the assertion.
Notably, the bound on the L 2 -norms in Theorem 7.5 is almost independent of the discretization parameters. Only the quadrature weights of the complexity reduction step enter. To make the bound uniform, one has to require the quadrature weights to be bounded from below by a positive constant. The concluding result now reads as follows.
Proof. Given the solution (ρ, m) exists up t, t > 0, we define F(t) := max i∈I ξ i − 1 2 [C 1 H c (ρ(t), m(t)) + C 2 ] with constants C 1 , C 2 as in Theorem 7.5. From the energy dissipation (Theorem 3.4) and the boundedness of the trace operator T , it follows d dt F(t) = C 3 d dt H c (ā(t)) ≤ C 3 u(t) · T m(t) ≤ C 4 ||u(t)|| ||m(t)|| for constants C 3 , C 4 > 0. Further, ||ρ(t)|| + ||m(t)|| ≤ F(t) holds by Theorem 7.5. Thus, it holds by the former estimate that d/dtF(t) ≤ C 4 ||u(t)||F(t). We follow by the Gronwall Lemma that F grows at most exponentially in time, and accordingly the same holds for the norm of the solution (ρ, m). By that existence of a continuously differentiable solution for any t > 0 can be deduced by the Peano existence theorem and the extension theorem for differential equations by standard arguments for ordinary differential equations, cf. [34].

Conclusion
We proposed a snapshot-based model reduction approach for a nonlinear flow problem on networks that is governed by the barotropic Euler equations with friction. The approach consists of a structure-preserving modification of the proper orthogonal decomposition combined with a quadrature-type complexity reduction. It yields online-efficient reduced models with remarkable structure-preserving properties, i.e., local mass conservation, an energy bound and a port-Hamiltonian structure. These properties are established using a few compatibility conditions, which we assure by appropriate adaptions in the training phase. The involved implications and appropriate efficient algorithmic implementations were a main focus of this paper. Notably, our model order reduction spaces fulfill an optimality condition under the compatibility conditions. As we demonstrated on a realistic gas transportation network benchmark, our reduced models show superior stability and overall performance compared to more generic non-structure preserving model reduction approaches. Moreover, a well-posedness result for them was established under a few additional assumptions.