Index-Aware Model-Order Reduction for a Special Class of Nonlinear Differential-Algebraic Equations

We extend the index-aware model-order reduction method to systems of nonlinear differential-algebraic equations with a special nonlinear term f(Ex),\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathbf{f}(\mathbf{E}\mathbf{x}),$$\end{document} where E\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathbf{E}$$\end{document} is a singular matrix. Such nonlinear differential-algebraic equations arise, for example, in the spatial discretization of the gas flow in pipeline networks. In practice, mathematical models of real-life processes pose challenges when used in numerical simulations, due to complexity and system size. Model-order reduction aims to eliminate this problem by generating reduced-order models that have lower computational cost to simulate, yet accurately represent the original large-scale system behavior. However, direct reduction and simulation of nonlinear differential-algebraic equations is difficult due to hidden constraints which affect the choice of numerical integration methods and model-order reduction techniques. We propose an extension of index-aware model-order reduction methods to a special class of nonlinear differential-algebraic equations without any kind of linearization. The proposed model-order reduction approach involves automatic decoupling of nonlinear differential-algebraic equations into nonlinear ordinary differential equations and algebraic equations, based on the decoupling of the linear differential equations obtained by ignoring the nonlinear term, thanks to an additional structural condition. This allows applying standard model-order reduction techniques to both parts without worrying about the index. The same procedure can also be used to simulate nonlinear differential-algebraic equations by standard integration schemes.


Introduction
We consider nonlinear differential-algebraic equations (DAEs) of the form; where f(Ex) ∈ R n and E is a singular matrix, A ∈ R n×n , B ∈ R n×m , C ∈ R ×n . The symbol denote time differentiation. x ∈ R n and y ∈ R are the state and output vectors, respectively. The input function u ∈ R m must be smooth enough, with the smoothness requirements depending on the index of the DAE. DAEs are known to be difficult to simulate and the level of difficulty is measured using index concepts such as differential index, tractability index, etc. The higher the index, the more difficult to simulate the DAE. Moreover, in practice, often such descriptor systems have very large n, compared to the number m of inputs and the number of outputs, which are typically small. Despite the ever increasing computational power, dynamic simulation using the system (1) is costly, see [1][2][3]. We are interested in a fast and stable prediction of the dynamics of DAE models, and therefore the application of model-order reduction (MOR) is vital. MOR aims to reduce the computational burden by generating reduced-order models (ROMs) that have lower computational cost to simulate, yet accurately represent the original large-scale system behavior. MOR replaces (1) by a ROM where E r , A r ∈ R r ×r , f r (E r x r ) ∈ R r , B r ∈ R r ×m and y r ∈ R , C r ∈ R ×r , such that the reduced-order of the state vector x r ∈ R r is r n. A good ROM should have small approximation error y − y r in a suitable norm . for a desired range of inputs u. There are many MOR methods for nonlinear systems such as proper orthogonal decomposition (POD), POD in conjunction with the discrete empirical interpolation method (POD-DEIM), see [3][4][5]. However, applying these MOR methods directly to DAEs leads to ROMs which are inaccurate or very difficult to simulate and sometimes have no solution, see [1,6]. It is a common practice to first convert nonlinear DAEs to ordinary differential equations (ODEs) by using index reduction (reformulation) techniques in order to be able to apply standard MOR methods for nonlinear systems such as POD. However, this index reduction may lead to drift-off effects or instabilities in the numerical solutions and may also depend on the structure of the nonlinear DAE. In [7], index-aware MOR (IMOR) methods were proposed to eliminate the index problem to allow employing standard techniques with ease. However, these methods were dedicated to linear DAEs. We propose an IMOR method for nonlinear DAEs of the form (1) which does not involve any kind of linearization. The starting point is a generalization of the index concept to this special class of nonlinear DAEs. The proposed generalization relies heavily on its linear counterpart together with some structural assumptions on the nonlinear term, and is consistent with the accepted definition of tractability index of nonlinear DAEs (see for instance [8]). The proposed IMOR approach is realized in two steps. The first step involves automatically decoupling the nonlinear DAEs into nonlinear differential and algebraic parts. Then, each part can be reduced separately using standard MOR techniques. The decoupled system generated from the first step can be used for numerical simulations by applying numerical integration on the ODE part and then solving the algebraic part.
The IMOR procedure can be applied to nonlinear DAEs of the form (1) with arbitrary index. Here we concentrate on index-1 systems, which can be used for a first assessment of the procedure. Nevertheless, we observe that index-1 DAEs can be found in relevant applications, such as in the description of gas transport networks.
The paper is organized as follows. In Sect. 2, we discuss the background of decoupling of DAEs and the tractability index. In Sect. 3, we propose the automatic decoupling of nonlinear DAEs of the form (1) using special projectors. In Sect. 4, we discuss the proposed IMOR method for nonlinear DAEs. In Sect. 5, we apply the proposed IMOR method to the nonlinear DAEs arising from gas transport networks. In the final section, we present some numerical examples illustrating the performance of the proposed method.

Decoupling of Linear Constant Coefficient DAEs
In this section, we repeat the procedure of decoupling linear DAEs and the theory it is based on, as this is the basis for the nonlinear decoupling.

Weierstraß Canonical Form
Our decoupling strategy was initially used to understand the underlying structure of linear constant coefficient DAEs via the Weierstraß canonical form [8,9]. Assuming f(Ex) = 0 and that the matrix pencil (E, A) is regular, (1) can be written as a Weierstraß-Kronecker canonical form which leads to an equivalent decoupled system where J ∈ R k×k and N ∈ R (n−k)×(n−k) is a nilpotent matrix with index μ. The vector is the i-th derivative of the input data. The control input matrices arẽ B 1 ∈ R k×m andB 2 ∈ R (n−k)×m . Subsystems (3a) and (3b) represent the inherited ODE and algebraic part, respectively, and the solutions of (1) can be obtained using are commonly known as the slow and fast parts of the solution, respectively. In [9], only complex-valued matrix pencils are considered and even for a real-valued pencil the matrix J in (3a) will be complex-valued in general, since in the Weierstraß canonical form J contains the (possibly complex) eigenvalues of the ODE part of the pencil. This can be circumvented by using the quasi-Weierstraß form introduced in [10]. This form guarantees that J and N have entries in the same field as those of (E, A). An index concept was introduced to classify different types of DAEs with respect to the difficulty arising in the theoretical and numerical treatment of a given DAE. "Index" is a notion used in the theory of DAEs for measuring the distance from a DAE to its related ODE. There are several definitions of a DAE index. The index μ in (3) is known as the Kronecker (nilpotency, differentiability) index. An equivalent decoupled system (3) shows the dependence of the solution of a linear DAE on the derivatives of the input function. In (3b), we can observe that the input function has to be at least μ − 1 times differentiable. The higher the index the more differentiations of the input data are involved. Since numerical differentiation is an unstable process, the index μ is a measure of numerical difficulty when solving the DAE. We can also observe that the initial conditionx 1 (0) of the differential part can be chosen arbitrary whilex 2 (0) has to satisfy hidden constraint Thus, DAE (1) has a unique classical solutions if x(0) = x 0 is consistent. The index problem affects the choice of numerical integration schemes strongly if standard numerical integration schemes are applied to DAEs directly without decoupling. This lead to the development of numerical integration schemes which were specifically designed for DAEs, see [8,11,12]. Hence, a promising way to solve and apply MOR to DAEs is to first split them into differential and algebraic parts, see [6]. According to [11], transforming linear DAEs into a Kronecker canonical form is numerically infeasible and it is restricted to linear DAEs. Due to this drawback other index concepts, such as the tractability index, differentiation index, etc, see [13], were proposed with each of them stressing different aspects of the DAE. It is worth mentioning that there are other techniques, such as the staircase form [14], which avoids the computation of eigenvalues and provides a numerically stable algorithm for the decoupling of matrix pencils. We also mention the quasi-Weierstraß form from [10], which makes use of so called Wong sequences and does not require the computation of eigenvalues.

Tractability Index
In this paper, we consider the tractability index introduced in [15] and its generalization in [16] defined as in Definition 1. For a thorough discussion of this concept we refer to [8].
Definition 1 (Tractability index) Given a regular matrix pair (E, A), we can define a matrix and projector chain by setting E 0 := E and A 0 := A, and where Q j are projectors onto Ker E j and P j = I − Q j . Then, if there exists an index γ such that E j is singular for all 0 ≤ j ≤ γ − 1, while E γ is non-singular, we say that γ is the tractability index of (E, A).
It is possible to prove that such an index γ exists for any regular matrix pair [8]. This index criterion does not depend on the special choice of the projector functions Q j , see [17]. The tractability index has gained a lot of attention since it can be calculated without the use of derivative arrays [18]. Hence, it is numerically feasible to compute the tractability index compared to computing the Kronecker index. This is the main tool that we chose in the decoupling of linear DAEs into their differential and algebraic parts, since it allows automatic decoupling procedures, see [7,13,16]. In order to decouple linear DAEs with index higher than one, so-called canonical projectors were introduced in [19] with additional constraints Based on these projectors, special projector bases were introduced leading to a decoupled system of the same dimension as the original, see [7]. A key step in forming the projectors in (4) is to find the initial projectors Q j spanning the nullspaces of the usually sparse E j . Standard ways of identifying the nullspace include the singular value decomposition (SVD) or alike, which does not utilize matrix patterns and can be expensive for large-size matrices.
One feasible way is to employ the sparse LU decomposition-based routine, called LUQ, see [20]. This routine was also used to construct the projector bases efficiently. This motivated the development of the index-aware MOR methods in [7]. In [7], equivalent explicit and implicit decoupling methods were proposed which are discussed in the following. According to [7], if we also assume f(Ex) = 0 and that the matrix pencil (E, A) is regular, using the special projectors proposed in [19] and projector bases introduced in [21], DAE (1) can be rewritten into an equivalent explicit decoupled system given by where L ∈ R n q ×n q is a nilpotent matrix with index γ. u ( j) and ξ ( j) p are the j-th derivatives with respect to t. The subsystems (6a) and (6b) correspond to the differential and algebraic parts of system (1). ξ p ∈ R n p and ξ q ∈ R n q are the differential and algebraic variables. The dimension of the decoupled system is given by n = n p + n q . The original variable x can be reconstructed by where M p ∈ R n×n p , M q ∈ R n×n q , and the overall matrix (M p , M q ) is invertible. Both M p and M q are constructed iteratively using a finite number of steps, and depend on the projector chain which leads to the index definition. We can observe that decoupled systems (3) and (6) are equivalent if A q = 0. Decoupled system (6) can be constructed automatically, and thus is numerically feasible. However, according to [7] the decoupling procedure of (6) involves the inversion of the non-singular matrix E γ which is costly for large-scale systems. Then, the implicit version of (6) was also proposed in [7] which does not involve the inversion of non-singular matrix E γ . Using this decoupling procedure, DAE (1) can be re-written into an equivalent implicit decoupled system given by where N q = L L −1 q is also a nilpotent matrix with the same index γ as L . The matrices L q ∈ R n q ×n q and E p ∈ R n p ×n p are always non-singular, see [21]. The subsystems (8a) and (8b) correspond to the differential and algebraic parts of system (1). ξ p ∈ R n p and ξ q ∈ R n q are the differential and algebraic variables. We can observe that the inherited ODEs (6a) and (8a) of the explicit and implicit decoupled systems can be simulated using standard ODE integration schemes. After obtaining the solutions of (8a), the algebraic part (8b) can be solved using numerical solvers such as LU decomposition-based routines. It is not straightforward to extend this to nonlinear DAEs. However, specific classes of nonlinear DAEs have been studied, usually those appearing in practice, namely DAEs with quadratic nonlinearity [22].

Decoupling of Nonlinear DAEs
In this section, we propose the decoupling of a class of nonlinear DAEs of the form (1a). This decoupling strategy is an extension of the decoupling strategy for linear DAEs proposed in [19]. A more general extension to nonlinear DAEs can be found in [8]. Here we adopt a simplified approach, which is feasible thanks to the special structure of the DAEs we are considering.

Decoupling Using Projectors
The tractability index of (1a) is independent of the nonlinearity if we assume the condition This implies that, after we apply the decoupling procedure for the linear part, the nonlinear term depends only on the differential variables. In [7], it was shown that (6) can be rewritten in the compact form Then using (11), we arrive at the explicit nonlinear decoupled system, The linear coefficients of the above system can be constructed in following the procedure in [7]. The nonlinear functions are defined as and The matrix coefficients A p f ∈ R n p ×n and A q f ∈ R n q ×n can also be constructed iteratively from projector bases proposed in [7]. This procedure can be applied to systems of arbitrary index. Here, we concentrate on index-1 systems, showing the explicit form of the additional condition (9).
can be written as We choose a projector Q 0 such that ImQ 0 = KerE 0 and its complementary projector P 0 = I − Q 0 . Using (4), which satisfy the identities: Substituting the above identities into (15) and simplifying leads to If we assume E 1 to be nonsingular, then (17) can be written as Since E 1 is nonsingular, then we say that the nonlinear DAE (1) is of tractability index 1. Left multiplying (18) by projectors P 0 and Q 0 separately, we obtain the differential and algebraic subsystems, respectively, of (1) given by where x P = P 0 x and x Q = Q 0 x. It is immediate to see that condition (9) is equivalent to EQ 0 = 0 in the index-1 case, which holds by construction. This procedure can be extended to higher index if projectors, constructed with the additional condition (5), are used. Then the original variable x can be decomposed as [7] x = x P + x Q , In order to ensure that the nonlinear term in (1a) depends only on the differential variable, that is, we need to enforce the following condition: which is equivalent to condition (9). A detailed discussion of condition (20) can be found in [23].
We can see that the decoupled system (19) is of dimension 2n while the DAE (1) is of dimension n. This implies that decoupling using projectors does not preserve the dimension of the original DAE. In the next section, we discuss how to derive a decoupled system which preserves the dimension of the nonlinear DAE (1).

Explicit Decoupling Using Bases
Projector bases can be applied to (19) as follows.
Let n q = dim(Ker E 0 ) and n p = n − n q . If, we also let q 0 ∈ Im Q 0 and p 0 ∈ Im P 0 , then, we can expand x with respect to the bases, obtaining where ξ q ∈ R n q , ξ p ∈ R n p , which implies that x P = p 0 ξ p and x Q = q 0 ξ q in (19). The left inverses of matrices q 0 ∈ R n×n q and p 0 ∈ R n×n p are denoted by q * T 0 ∈ R n q ×n and p * T 0 ∈ R n p ×n , respectively. Substituting x P = p 0 ξ p and x Q = q 0 ξ q into (19) leads to a decoupled system which can be left multiplied by the left inverses p * T 0 and q * T 0 , respectively. This yields a decoupled system in compact form: where Using the modified decoupling procedure, we can identify the matrices in (7) with M p = p 0 , M q = q 0 , so that condition (9) is equivalent to Eq 0 = 0, which is satisfied by construction. We can observe that system (22) has the same structure as system (12) with matrix coefficients given by (10).
We can now observe that the total dimension of the decoupled system is n = n p + n q , which is equal to the dimension of the nonlinear DAE (1). Instead of solving the coupled nonlinear DAE (1) we can now solve the decoupled nonlinear system (22). We obtain the solution ξ p by applying standard integration schemes to (22a) and the solutions of ξ q can be computed by post-processing using (22b). Then, the desired output solution can be obtained using (22c). However, we can observe that the coefficients of (22) involve computing the inverse of E 1 which is computationally expensive and requires large storage for large scale systems. Moreover, it also leads to dense matrix coefficients of the decoupled system (22).

Implicit Decoupling
In this subsection, we discuss a decoupling strategy which does not involve inversion of matrix E 1 . In [7], it was shown that (8) can be rewritten in the compact form where L q and L are as defined in (8). Then using the above system we can deduce the implicit nonlinear decoupled system of (1) with condition (9) given by where nonlinear functions f p (ξ p ) ∈ R n p and f q (ξ p ) ∈ R n q are defined as in (13) and (14). The construction of the matrix coefficients of system (24) can be done using the same procedure as for the case of linear DAEs in [7] for index greater than one. For simplicity, we discuss the index-1 systems. This is done as follows. Substituting (21) into (17) leads to Instead of inverting matrix E 1 , we can decouple (25) into differential and algebraic parts using matricesp 0 ∈ R n×n p andq 0 ∈ R n×n q proposed in [7] which are defined as viâ The system (26) can be reduced to a nonlinear decoupled system given by where The nonlinear terms are defined as: We note that matrices E p and E q are always nonsingular. As expected, system (27) is equivalent to system (24) with matrix coefficients given by It is also clear that f(E 1 p 0 ξ p ) = f(Ep 0 ξ p ), in accordance with (10). We can observe that (27) does not involve any matrix inversions. It is an implicit version of the decoupled system (22) and their output solutions must coincide. However, in practice it is computationally cheaper to construct the coefficients of (27) than those in (22). Both decoupled systems preserve the dimension and the stability of the nonlinear DAE (1). If (1) is of tractability index 1, then it can be automatically decoupled into either (27) or (22). Thus, instead of simulating (1), we can simulate its equivalent nonlinear decoupled system (27) easily using standard numerical integration and solvers. Decoupled systems (22) and (27) can be constructed in efficient way by employing the sparse LU decomposition-based routine, called LUQ, see [20], to construct the projectors and their respective bases. In the next section, we discuss how to apply MOR to (27).

Index-Aware MOR for Nonlinear DAEs
Here, we consider the equivalent nonlinear decoupled system (24) corresponding to the nonlinear DAE (1), but the same strategy can be applied to (12). Given such a nonlinear decoupled system, our goal is to reduce the order of differential and algebraic parts separately.

MOR for the Nonlinear Differential Subsystem
We consider the nonlinear differential subsystem of the nonlinear decoupled system (24) given by where y p ∈ R is the output solution of the differential part. Our goal is reduction by projection of system (28). This means we want to find a linear subspace in which the solution trajectory lies approximately. This subspace is defined by its basis matrix V p ∈ R n p ×r p where r p n p . We are interested in finding a solution ξ p r ∈ R r p such that ξ p ≈ V p ξ p r . We can then project system (28) onto that subspace by Galerkin projection resulting in the reduced differential subsystem where Projection matrix V p can be computed using standard MOR techniques for nonlinear systems such as POD [24]. However, if we employ model order reduction via POD by using (29a) to compute the snapshots, the nonlinearity f p r (ξ p r ) = V T p f p (V p ξ p r ) requires computation of f p (V p ξ p r ) which has a complexity in the system dimension. Therefore, the discrete empirical interpolation method (DEIM) [5,25] is often used to create a truly low-dimensional function approximating V T p f p (V p ξ p r ). In this paper we do not consider DEIM or other hyperreduction methods.

Reduction of Algebraic Subsystem
After reducing the differential subsystem using, for example, POD, the nonlinear term in the algebraic subsystem (24b) is also affected leading to where y q ∈ R is the output solution of the algebraic part after reducing the differential subsystem. Here, we intend to reduce the size of the algebraic variables ξ q by constructing another matrix V q ∈ R n q ×r q where r q n q . That is, we replace (30) by a reduced algebraic subsystem given by −L r ξ q r =A q r ξ p r −L q r ξ q r + f q r (ξ p r ) + B q r u, (31a) Reduction matrix V q can also be computed using the POD by taking the algebraic solutions of (24b) obtained by solving the linear system where N q is as defined in (8).
Combining (29) and (31), we obtain an index-aware reduced order model (I-ROM) of (1) given by where the reduced dimension is given by r = r p + r q n. Thus, we replace (1) with (33) instead of (2).

Nonlinear DAEs Arising from Gas Networks
In this section, we apply the implicit decoupling strategy proposed in Sect. 3.3 to nonlinear DAEs arising from gas flow in pipeline networks.

Index Reduction of DAEs Arising from Gas Networks
We consider a spatial discretization approach of one dimensional isothermal Euler equations arising from gas flow pipe networks proposed in [1,26], leading to a nonlinear DAE given by The unknowns are described by the pressure at the supply nodes p s ∈ R n s , the pressure at all other nodes p d ∈ R n d +n 0 , the difference of flux over a pipe segment q − ∈ R n E and the average of the mass flux over a pipe segment q + ∈ R n E , modelled over a graph with n E edge segments, that correspond to the size of the discretization, n s supply nodes, n d demand nodes and n 0 interior nodes. The diagonal matrices M L ∈ R n E ×n E and M A ∈ R n E ×n E encode parameters such as length, radius of the pipe segments as well as constants coming from the gas equation. The matrix B d ∈ R (n d +n 0 )×n d is a matrix of ones and zeros making sure that the demand of the demand node is put at the right place in the mass flux equation. The matrix A 0 ∈ R n d ×n E is extracted from the incidence matrix of the graph representing the refined gas transportation network and removing the rows corresponding to the supply nodes, while A S ∈ R n s ×n E is the matrix extracted from the incidence matrix by only taking rows corresponding to the supply nodes.The absolute value sign around the matrices |A 0 | and |A S | mean that component-wise absolute values are taken to form the matrix, see [1].  , s i (t), . . .) T ∈ R m s are vectors for flux (mass flow) at demand nodes and pressure at supply nodes, respectively. The nonlinear term is the vector involving friction and gravitation effects with where ψ k (p d , p s ) is the k-th entry of the vector-valued function: The scalars λ k , D k , L k and A k denote friction, diameter, length and area of the pipe's k-th segment. The scalar h k denotes the height difference of the pipe segment. These scalar parameters in the system and those defined earlier are known at least within some range of uncertainty. System (34) can be rewritten in the form (1) leading to a system of nonlinear DAEs with dimension n = 2n E + n d + n 0 + n s . The desired outputs in R n s +n d can be obtained using the output equation where y q = |A S |q + is the mass flow at the supply nodes and y p = B T d p d is the pressure at demand nodes. We can observe that the initial condition has to be consistent with the hidden constraints in (34). Efficient simulation of (34) has numerical integration challenges since the solutions of hyperbolic balance laws can blow-up in finite time, due to both the stiffness and index problem. In [26], an index reduction strategy was proposed to eliminate the index problem. This was done by reformulating (34) into an implicit nonlinear ODE given by .
Since from (36) we are just interested in the solutions of q + and p q , the dimension of the nonlinear DAE (34) can be reduced toñ = n d + n 0 + n E with output equation The generated ODE can be reduced further using standard MOR methods for nonlinear systems, such as POD, POD-DEIM, etc, applied to (37), see [1]. However, the index reduction approach presented depends on the spatial discretization approach used.
In the next section, we propose an alternative model which preserves the DAE structure independent of the spatial discretization method.

Decoupled Model of Gas Transport Networks
Here, we discuss the decoupling analysis of nonlinear DAE (34) arising from the gas transportation networks. We can observe that (34) can be re-written into the form (1) where The unknown vector x and input vector u are given by Since the gas transport model can be rewritten in the form (1), we can decoupled it into either the form (12) or (23). In our discussion, we shall use the implicit decoupling strategy proposed in Sect. 3.3 leading to an implicit decoupled system (23). For convenience, we can partition (38) into a block form leading to ⎛ where The nonlinear term is defined as (40) From the above expression, we can observe that the nonlinear function in (39a) only depends on the differential variables, by noting that the algebraic variables p s in x 3 are not included in the definition ofg. This implies that the tractability index of this system is independent of the nonlinearity which satisfies condition (9). Hence we can apply our proposed decoupling procedure in Sect. 3. In order to decouple (39), we need to first find the tractability index of (39) using Definition 1. Setting we can then construct projectors such that E 0 Q 0 = 0, meaning E 13 Q = 0, or Q ∈ R n v ×n v is the projector onto the nullspace of E 13 and P ∈ R n v ×n v is its complementary projector. Substituting the above matrices and projectors into (4) leads to If E 1 is nonsingular, the DAE (39) is of tractability index 1. Next, we construct the values of the matrix coefficients of (27) as follows. Let n p = rank(E 0 ) and n q = n − n p . Then, the columns of the matrices are linearly independent and span the column spaces of Q 0 and P 0 in (42), respectively. The left inverse of matrices p 0 and q 0 are given by respectively, where q * T and p * T are the left inverses of column matrices q and p, respectively. Let k q be the dimension of the nullspace of E 13 , and k p = n v −k q . The columns of q ∈ R n v ×k q and p ∈ R n v ×k p are linearly independent and span the column spaces of Q and P in (42), respectively. Finally, column matricesp 0 ∈ R n×n p andq 0 ∈ R n×n q can be constructed such that their columns are linearly independent and span the null spaces of the matrices q T 0 A T 0 ∈ R n q ×n and E T 0 ∈ R n×n , respectively. The differential and algebraic variables are given by 3 and ξ q = q * T 0 Q 0 x = respectively. The nonlinear term is defined as andẼ 13 = E 13 p. This is due to the fact that E 13 x 3 = E 13 pp * T x 3 = E 13 pξ p 2 . It can be proved that f q (ξ p ) =q T 0f (ξ p ) = 0 ∈ R n q always due to the structure of the nonlinearity. Finally, substituting (41), (43)-(45) into (27) leads to an equivalent nonlinear decoupled system of (38) given by where f p (ξ p ) =p T 0f (ξ p ) ∈ R n p withf(ξ p ) as defined in (46). The system matrix coefficients are computed as defined in (27). We can observe that the nonlinear gas transport network model has been decoupled into n p = n E + k p nonlinear differential equations, and n q = n E + k q algebraic equations. Subsystem (47a) can be simulated using standard numerical integration, then algebraic solutions of (47b) can be obtained by using numerical solvers after post-processing. Hence, the desired output data can be obtained through (47c). We note that the decoupling enables us to treat DAEs like ODEs. However, the stiffness problem is inherited in the ODE subsystem (47a). In order to cope with the stiffness problem, we can use IMEX integration scheme [27] instead of standard integration which makes an efficient simulation of (47a) possible. We note that the values of the matrix coefficients of (47) can vary depending on the choices of projectors in (42), but the solutions will always be the same. In practice, system (47) can be constructed automatically following the implicit decoupling procedure in Sect. 3.3. Numerical experiments show that (47a) and (37) have the same dimension for the case of index 1 gas transportation networks.

Numerical Experiments
In this section, we illustrate the performance of the proposed decoupling and IMOR method for nonlinear DAEs with a special nonlinear term f(x) = f(Ex), where E is a singular matrix. Such nonlinear DAEs can arise from gas transportation networks as discussed in Sect. 5. Here, we consider small to large examples of gas transportation networks leading to nonlinear DAEs of tractability index 1. We compute the relative error in the format Re.error = y−y r 2 / y 2 . The output error is defined as max(Re.error(pressure), Re.error(mass flow)).
Simulations were done using MATLAB®Version 2012b on a Unix desktop.

Numerical Integration
We compare the output solutions (mass flow at the supply node and pressure at demand nodes) of different gas transportation models: nonlinear DAE model (34), nonlinear ODE model (37) and nonlinear decoupled model (47).

Example 1
In this example, we consider small to medium gas pipeline networks from [28,29] with steady pressure at the supply pressure node and steady mass flow at demand nodes. We are interested in the comparison of the pressure and mass flows of different models of each gas transportation network shown in Table 1.
In Table 1, we can observe that the index reduced ODE model has the same dimension as the differential part of the nonlinear decoupled model. We use the implicit-Euler numerical integration scheme to solve the nonlinear DAE and ODE models with a fixed time step. For the nonlinear decoupled model we use the implicit-Euler numerical integration scheme on the differential part and LU based numerical solver for the algebraic part. Figures 1, 2, 3, 4, show the mass flow at the supply node and pressure at the first demand node for each network presented in Table 1. In Fig. 1, we used steady pressure s(t) = 650bars at the supply node and steady mass flow rate of d(t) = 100Kg/s at the demand node.
In Fig. 2, we used steady pressure s(t) = 700bars at the supply node and steady mass flow rate of d(t) = (60, 30) T at the demand nodes.
In Fig. 3, we used steady pressure s(t) = 4.55 × 10 3 bars at the supply node and steady mass flow rate of d(t) = (0. 21 at the demand nodes. In Fig. 4, we used steady pressure s(t) = 3.45 × 10 4 bars at the supply node and steady mass flow rate of d(t) = 10 × ones(24, 1) at the demand nodes. In all test cases, we can observe that all models decay towards steady mass flow at the supply node and steady pressure at the demand nodes.

Example 2
In this example, we are interested in comparing the pressure and mass flow rate while applying steady pressure at supply node and transient mass flow rate at the demand node. We consider a medium size gas transport pipe network with 200 pipes, one supply node and one demand node generated using the following data. The length, diameter and average roughness of each pipe are chosen as constants given by 18.15m, 1.422m and 1.5 × 10 −6 m, respectively. The gas composition through the network is methane with specific gas constant 518.26J/KgK at steady supply of 84bar and mass flow at demand as shown in the first row of Fig. 5 in the time interval t ∈ [0, 1000s] .
This leads to a nonlinear DAE system of dimension n = 601 which we decoupled into n p = 400 differential equations and n q = 201 algebraic equations. For comparison, we    In this example, we compare the matrix properties of the matrix pencils of the derived models and the values of the nonlinear term at a fixed state vector. In Figs. 6, 7, 8, we compare the sparsity of the matrix pencils of the coupled model, decoupled model and implicit ODE model. We can observe that all models are sparse, however the decoupled model is the least sparse. In Table 2, we compare the finite spectrum of the matrix pencils and the nonlinearity. For the definition of finite spectrum and singular values of a matrix pencil, see [8]. We can observe that all models have the same spectrum with purely imaginary finite eigenvalues and approximately the same values of the nonlinear function.

Model Order Reduction
Here, we illustrate the performance of the proposed IMOR method compared to existing MOR methods.

Example 4
We consider a large-scale gas transport pipeline network with 5,000 pipes, 1 supply node and 1 demand node. This model was generated numerically using the following data. The length, diameter and average roughness of each pipe are chosen as 0.726m, 1.422m and 1.0×10 −6 m, respectively. The gas composition is with specific gas constant 1530J/KgK at steady pressure 50bar at supply node and mass flow as a step function as shown in the first row of Fig. 10 at the demand node at a time interval t ∈ [0, 86400] . This leads to a nonlinear DAE (34) of dimension n = 15, 001. It took 63.7s to automatically decouple the nonlinear DAE into n p = 10, 000 nonlinear differential equations and n q = 5, 001 algebraic   equations. We also generated an index reduced ODE (37) of dimensionñ = 10, 000. We reduced the decoupled system using POD on both the differential and algebraic parts.
Then, we obtained an I-POD model with r p = 2 and r q = 4 leading to a total reduction of r = r p + r q = 6 15, 001. The first singular values of the differential part are shown in Fig. 9. We also used POD to reduce both the nonlinear DAE and ODE directly. For comparison, the size of ROMs for different MOR methods is determined by making sure that the output error is below 10 −4 and the results are presented in Table 3. All numerical integration was done using the backwards Euler method with a fixed time step h = 250 and the LU-based numerical solver implemented in the backslash MATLAB command was used for solving the resulting linear systems of equations in each time step. We note that using the backslash operator with a sparse matrix [30] in MATLAB, the LU-decomposition is done by UMF-PACK, i.e., row and column permutations are not only applied for numerical stability but also to reduce the number of fill-ins significantly. We can observe that I-POD leads to the largest ROM and lowest speed-ups. This is due to the fact that its ROM is a DAE while the other ROMs are ODEs. The comparison of the mass flow at the supply node and the pressure at the demand node of all ROMs are shown in Figs. 10 and 11, where the training and the test scenarios are different in Fig. 11. In Fig. 12, we compare the output relative error for pressure 6 8  Fig. 10. We can observe that I-POD is the most accurate while ODE-POD is the least accurate. However, I-POD leads to a slightly bigger ROM.

Conclusions
We have proposed a new automatically decoupling strategy and an IMOR method for nonlinear DAEs with a special nonlinear term. This approach eliminates the index problem during simulation and MOR which allows the use of standard numerical integration methods and MOR techniques. We have derived both the implicit (12) and explicit (23) decoupled systems for arbitrary index nonlinear DAEs. We have demonstrated the accuracy of this approach by applying it to nonlinear DAEs arising from the gas transportation networks. The computational cost of this approach can be improved by applying reordering algorithms after decoupling. However, we have restricted ourselves to nonlinear DAEs arising from gas transportation networks. Future research will deal with nonlinear DAEs arising from other applications. Funding Open Access funding enabled and organized by Projekt DEAL.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.