# Automatic Decoupling and Index-aware Model-Order Reduction for 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), where E 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 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. 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 using standard integration schemes. We illustrate the performance of our proposed method for nonlinear differential-algebraic equations arising from gas flow models in pipeline networks.

## Authors

• 1 publication
• 1 publication
• 7 publications
• 57 publications
07/16/2021

11/18/2020

### Continuous Galerkin Schemes for Semi-Explicit Differential-Algebraic Equations

This paper studies a new class of integration schemes for the numerical ...
01/05/2022

### An efficient extended block Arnoldi algorithm for feedback stabilization of incompressible Navier-Stokes flow problems

Navier-Stokes equations are well known in modelling of an incompressible...
11/24/2020

### Model Order Reduction for Gas and Energy Networks

To counter the volatile nature of renewable energy sources, gas networks...
07/10/2019

### Improved Structural Methods for Nonlinear Differential-Algebraic Equations via Combinatorial Relaxation

Differential-algebraic equations (DAEs) are widely used for modeling of ...
11/25/2020

### Towards a reliable implementation of least-squares collocation for higher-index differential-algebraic equations

In this note we discuss several questions concerning the implementation ...
12/25/2021

### Simulation of crumpled sheets via alternating quasistatic and dynamic representations

In this work, we present a method for simulating the large-scale deforma...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

We consider nonlinear differential-algebraic equations (DAEs) of the form;

 Ex′ =Ax+f(Ex)+Bu,Ex(0)=Ex0, (1a) y =Cx, (1b)

where and is a singular matrix, . The symbol denote time differentiation. and

are the state and output vectors, respectively. The input function

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 compared to the number 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 [10, 2, 6]. 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

 Erx′r (2a) yr =Crxr, (2b)

where and such that the reduced-order of the state vector is A good ROM should have small approximation error in a suitable norm for a desired range of inputs

There exist many MOR methods for nonlinear systems such as proper orthogonal decomposition (POD), POD in conjunction with the discrete empirical interpolation method (POD-DEIM), see

[6]. 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 [10, 4]. 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 [3], 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 index-aware MOR (IMOR) method for nonlinear DAEs of the form (1) which does not involve any kind of linearization. This 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 paper is organized as follows. In Section 2, we discuss the background of decoupling of DAEs and the tractability index. In Section 3, we propose the automatic decoupling of nonlinear DAEs of the form (1) using special projectors. In Section 4, we discuss the proposed IMOR method for nonlinear DAEs. In Section 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.

## 2 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.

### 2.1 Weierstraß canonical form

Our decoupling strategy was initially used to understand the underlying structure of linear constant coefficient DAEs via the Weierstraß canonical form [13]. Assuming and that the matrix pencil is regular, (1) can be written as a Weierstraß-Kronecker canonical form which leads to an equivalent decoupled system

 ~x′1 =J~x1+~B1u,~x1(0)=~x10, (3a) ~x2 =−μ−1∑i=0Ni~B2u(i), (3b)

where and is a nilpotent matrix with index The vector is the -th derivative of the input data. The control input matrices are and Subsystems (3a) and (3b) represent the inherited ODE and algebraic part, respectively , and the solutions of (1) can be obtained using where is a nonsingular matrix and The vectors and

are commonly known as the slow and fast parts of the solution, respectively. 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 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 condition of the differential part can be chosen arbitrary while has to satisfy hidden constraint

 ~x2(0)=−μ−1∑i=0Ni~B2u(i)(0).

Thus, DAE (1) has a unique classical solutions if 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 [15, 14]. Hence, a promising way to solve and apply MOR to DAEs is to first split them into differential and algebraic parts, see [4]. According to [15], 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 [20], were proposed with each of them stressing different aspects of the DAE.

### 2.2 Tractability index

In this paper, we consider the tractability index introduced in [8] and its generalization in [17] defined as in Definition 2.1.

###### Definition 2.1 (Tractability index ([17]) ).

Given a regular matrix pair We define a matrix and projector chain by setting and , given by

 Ej+1 :=Ej−AjQj,Aj+1:=AjPj,forj≥0, (4)

where are projectors onto and . There exists an index such that is non-singular and all are singular for all . is the tractability index of a DAE.

This index criterion does not depend on the special choice of the projector functions , see [18]. The tractability index has gained a lot of attention since it can be calculated without the use of derivative arrays [19]. Hence, it is numerically feasible to compute the tractability index compared to computing the Kronecker index. This is the main tool in the decoupling of linear DAEs into their differential and algebraic parts, since it allows automatic decoupling procedures, see [17, 3, 20]. In order to decouple linear DAEs with index higher than one , so -called canonical projectors were introduced in [16] with additional constraint s , for Based on these projectors, special projector bases were introduced leading to a decoupled system of the same dimension as the original. A key step in forming the projectors in (4) is to find the initial projectors spanning the nullspaces of the usually sparse

. 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

[22]. This routine was also used to construct the projector bases efficiently. This motivated the development of the index-aware MOR methods in [3]. In [3], equivalent explicit and implicit decoupling methods were proposed which are discussed in the following. According to [3], if we also assume and that the matrix pencil is regular, using the special projectors proposed in [16] and projector bases introduced in [1], DAE (1) can be rewritten into an equivalent explicit decoupled system given by

 ξ′p =Apξp+Bpu,ξp(0)=ξp0, (5a) ξq =γ−1∑j=0Lj(Aqξ(j)p+Bqu(j)) (5b) y =Cpξp+Cqξq, (5c)

where is a nilpotent matrix with index and are the -th derivatives with respect to . The subsystems (5a) and (5b) correspond to the differential and algebraic parts of system (1). and are the differential and algebraic variables. The dimension of the decoupled system is given by We can observe that decoupled systems (3) and (5) are equivalent if Decoupled system (5) can be constructed automatically, and thus is numerically feasible. However, according to [3] the decoupling procedure of (5) involves the inversion of non-singular matrix which is costly for large-scale systems. Then, the implicit version of (5) was also proposed in [3] which does not involve the inversion of non-singular matrix . Using this decoupling procedure, DAE (1) can be re-written into an equivalent implicit decoupled system given by

 Epξ′p =Apξp+Bpu,ξp(0)=ξp0, (6a) Lqξq =γ−1∑j=0Njq(Aqξ(j)p+Bqu(j)), (6b) y =Cpξp+Cqξq, (6c)

where is also a nilpotent matrix with the same index as The matrices and are always non-singular , see [1]. The subsystems (6a) and (6b) correspond to the differential and algebraic parts of system (1). and are the differential and algebraic variables. We can observe that the inherited ODEs (5a) and (6a) of the explicit and implicit decoupled systems can be simulated using standard ODE integration schemes. After obtaining the solutions of (6a), the algebraic part (6b) can be solved using numerical solvers such as LU decomposition-based routines. It is not straight forward to extend this to nonlinear DAEs. However, specific classes of nonlinear DAEs have been studied, usually those appearing in practice, see [5].

## 3 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 [16].

### 3.1 Decoupling using projectors

Assume that the tractability index of (1a) is independent of the nonlinearity, i.e., all projectors constructed using Definition 2.1 are constant matrices. Setting (1a) can be written as

 E0x′ =A0x+f(E0x)+Bu. (7)

We choose a projector such that and its complementary projector Using (4),

 E1=E0−A0Q0,A1=A0P0,

which satisfy the identities :

 E1P0 =E0,A1−E1Q0=A0. (8)

Substituting the above identities into (7) and simplifying leads to

 E1[P0x′+Q0x] =A1x+f(E1P0x)+Bu. (9)

If we assume to be nonsingular, then (9) can be written as

 P0x′+Q0x =E−11[A1x+f(E1P0x)+Bu]. (10)

Since is nonsingular, then we say that the nonlinear DAE (1) is of tractability index Left multiplying (10) by projectors and separately, we obtain the differential and algebraic subsystems, respectively, of (1) given by

 x′P =P0E−11A0xP+P0E−11f(E1xP)+P0E−11Bu,xP(0)=P0x(0), (11a) xQ =Q0E−11A0xP+Q0E−11f(E1xP)+Q0E−11Bu, (11b) y =CxP+CxQ, (11c)

where and We can see that decoupled system (11) is of dimension while the DAE (1) is of dimension . 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).

### 3.2 Explicit decoupling using bases

Projector bases can be applied to (11) as follows. Let and If, we also let and then, we can expand with respect to the bases, obtaining

 x=q0ξq+p0ξp, (12)

where which implies that in (11). The left inverses of column matrices and are denoted by and respectively. Substituting
into (11) leads to a decoupled system which can be left multipl ied by the left inverses and respectively. This yields a decoupled system in compact form:

 ξ′p =Apξp+fp(ξp)+Bpu,ξp(0)=p∗T0x(0), (13a) ξq =Aqξp+fq(ξp)+Bqu, (13b) y =Cpξp+Cqξq, (13c)

where

 Ap Cq =Cq0∈Rℓ×nq,Aq=q∗T0E−11A0p0∈Rnq×np,Bq=q∗T0E−11B∈Rnq×m.

and

 fp(ξp)=p∗T0E−11f(E1p0ξp)∈Rnp,fq(ξp)=q∗T0E−11f(E1p0ξp)∈Rnq.

We can now observe that the total dimension of the decoupled system is , 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 (13). We obtain the solution by applying standard integration schemes to (13a) and the solutions of can be computed by post-processing using (13b). Then, the desired output solution can be obtained using (13c). However, we can observe that the coefficients of (13) involve computing the inverse of which is computationally expensive and requires large storage for large scale systems. Moreover, it also leads to dense matrix coefficients of the decoupled system (13).

### 3.3 Implicit decoupling

In this subsection, we discuss a decoupling strategy which does not involve inversion of matrix . This is done as follows. Substituting (12) into (9) leads to

 =(A0p0−E1q0)(ξpξq)+f(E1p0ξp)+Bu. (14)

we can decouple (14) into differential and algebraic parts using column matrices and proposed in [3] which are defined as via and . Left multiplying (14) by leads to

 (^pT0E1p0000)(ξpξq)′=(^pT0A0p00^qT0A0p0−^qT0E1q0)(ξpξq)+(^pT0f(E1p0ξp)^qT0f(E1p0ξp))+(^pT0B^qT0B)u. (15)

The system (15) can be reduced to a nonlinear decoupled system given by

 Epξ′p =Apξp+fp(ξp)+Bpu,ξp(0)=p∗T0x(0), (16a) Eqξq =Aqξp+fq(ξp)+Bqu, (16b) y =Cpξp+Cqξq, (16c)

where

 Ep =^pT0E0p0∈Rnp×np,Ap=^pT0A0p0∈Rnp×np,Bp=^pT0B∈Rnp×m, Eq =−^qT0A0q0∈Rnq×nq,Aq=^qT0A0p0∈Rnp×nq,Bq=^qT0B∈Rnq×m.

The nonlinear terms are defined as: where We note that matrices and are always nonsingular. We can observe that (16) does not involve any matrix inversions. It is an implicit version of the decoupled system (13) and their output solutions must coincide. However, in practice it is computationally cheaper to construct the coefficients of (16) than those in (13). Both decoupled systems preserve the dimension and the stability of the nonlinear DAE (1). If (1) is of tractability index , then it can be automatically decoupled into either (16) or (13). Thus, instead of simulating (1), we can simulate its equivalent nonlinear decoupled system (16) easily using standard numerical integration and solvers. Decoupled systems (13) and (16) can be constructed in efficient way by employing the sparse LU decomposition-based routine, called LUQ, see [22], to construct the projectors and their respective bases. In the next section, we discuss how to apply MOR to (16).

## 4 Index-aware MOR for nonlinear DAEs

Here, we consider the equivalent nonlinear decoupled system (16) corresponding to the nonlinear DAE (1), but the same strategy can be applied to (13). Given such a nonlinear decoupled system, our goal is to reduce the order of differential and algebraic parts separately.

### 4.1 MOR for the nonlinear differential subsystem

We consider the nonlinear differential subsystem of the nonlinear decoupled system (16) given by

 Epξ′p =Apξp+fp(ξp)+Bpu,ξp(0)=p∗T0x(0), (17a) yp =Cpξp, (17b)

where is the output solution of the differential part. Our goal is reduction by projection of system (17). This means we want to find a linear subspace in which the solution trajectory lies approximately. This subspace is defined by its basis matrix where We are interested in finding a solution such that We can then project system (17) onto that subspace by Galerkin projection resulting in the reduced differential subsystem

 Eprξ′pr =Aprξpr+fpr(ξpr)+Bpru, (18a) ypr =Cprξpr, (18b)

where
and Projection matrix can be computed using standard MOR techniques for nonlinear systems such as POD [7]. However, if we employ POD by using (18a) to compute the snapshots, the nonlinearity requires computation of which has a complexity in the system dimension. Therefore, we use discrete empirical interpolation method (DEIM) to create a truly low-dimensional function approximating The DEIM algorithm creates matrices such that

 VTpfp(Vpξpr)≈VTpUp(WTpUp)−1WTpfp(Vpξpr).

Here is orthonormal and the matrix is a picking matrix, where each row has exactly one nonzero entry which is This means that picks functions from the vector of functions . Here we have to make sure to pick appropriately, in order to make truly low-dimensional, see [11].

### 4.2 Reduction of algebraic subsystem

After reducing the differential subsystem using, for example, POD, the nonlinear term in the algebraic subsystem (16b) is also affected leading to

 Eqξq ≈AqVpξpr+fq(Vpξpr)+Bqu, (19a) yq ≈Cqξq, (19b)

where is the output solution of the algebraic part after reducing the differential subsystem. Here, we intend to reduce the size of the algebraic variables by constructing another matrix where That is, we replace (19) by a reduced algebraic subsystem given by

 Eqrξqr =Aqrξp+fqr(ξpr)+Bqru, (20a) yqr =Cqrξqr, (20b)

where
and Reduction matrix can also be computed using the POD by taking the algebraic solutions of (16b) obtained from the snapshots of (16b) as snapshots. Also here, the nonlinearity has to be evaluated completely, even though we reduce the algebraic system size. Hence, we also need to use the DEIM to create a truly low-dimensional function approximating The DEIM algorithm creates matrices such that

 VTqfq(Vqξpr)≈VTqUq(WTqUq)−1WTqfq(Vpξpr),

where and is a picking matrix. Combining (18) and (20a), we obtain an index-aware reduced order model (I-ROM) of (1) given by

 Eprξ′pr =Aprξpr+fpr(ξpr)+Bpru,ξpr(0)=ξpr0, (21) Eqrξqr =Aqrξp+fqr(ξpr)+Bqru, yr =Cprξpr+Cqrξqr,

where the reduced dimension is given by Thus, we replace (1) with (21) instead of (2).

## 5 Nonlinear DAEs arising from gas networks

In this section, we apply the implicit decoupling strategy proposed in Subsection 3.3 to nonlinear DAEs arising from gas flow in pipeline networks.

### 5.1 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 [12, 10], leading to a nonlinear DAE given by

 =−M−1Lq−, (22a) ∂tq+ (22b) 0 =A0q++|A0|q−−Bdd(t), (22c) 0 =ps−s(t). (22d)

The unknowns are described by the pressure at the supply nodes , the pressure at all other nodes the difference of flux over a pipe segment and the average of the mass flux over a pipe segment , modelled over a graph with edge segments, that correspond to the size of the discretization, supply nodes, demand nodes and interior nodes. The diagonal matrices and encode parameters such as length, radius of the pipe segments as well as constants coming from the gas equation. The matrix 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 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 is the matrix extracted from the incidence matrix by only taking rows corresponding to the supply nodes. and are the incidence matrices of the undirected graph defined as the component-wise absolute values of the incidence matrices of the directed graph, see [10]. The input functions and 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

 gk(q+,pd,ps) =−gAk2γ0ψk(pd,ps)ΔhkLk−λkγ04DkAkqk+|qk+|ψk(pd,ps), (23)

where is the k-th entry of the vector-valued function:

 ψ(pd,ps)=|ATS|ps+|AT0|pd∈RnE.

The scalars and denote friction, diameter, length and area of the pipe’s -th segment. The scalar 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 (22) can be rewritten in the form (1) leading to a system of nonlinear DAEs with dimension The desired outputs in can be obtained using the output equation

 y =(yqyp)=(0|AS|0000BTd0)⎛⎜ ⎜ ⎜ ⎜⎝q−q+pdps⎞⎟ ⎟ ⎟ ⎟⎠, (24)

where is the mass flow at the supply nodes and is the pressure at demand nodes. We can observe that the initial condition has to be consistent with the hidden constraints in (22). Efficient simulation of (22) 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 [12], an index reduction strategy was proposed to eliminate the index problem. This was done by reformulating (22) into an implicit nonlinear ODE given by

 (|A0|ML|AT0|00I)(∂tpd∂tq+)=(0A0MAAT00)(pdq+)+(|A0|ML|ATS|∂ts(t)g(q+,s(t),pd))+(0−BdMAATS0)(s(t)d(t)). (25)

Since from (24) we are just interested in the solutions of and the dimension of the nonlinear DAE (22) can be reduced to with output equation

 y=(yqyp)=(0|AS|BTd0)(pdq+).

The generated ODE can be reduced further using standard MOR methods for nonlinear systems, such as POD, POD-DEIM, etc, applied to (25), see [10]. 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.

### 5.2 Decoupled model of gas transport networks

Here, we discuss the decoupling analysis of nonlinear DAE (22) arising from the gas transportation networks. As a result, we present an alternative model to the ODE model (25) proposed in [10]. We can observe that (22) can be re-written into the form (1) where

 E =⎛⎜ ⎜ ⎜ ⎜⎝00|AT0||ATS|0I0000000000⎞⎟ ⎟ ⎟ ⎟⎠,A=⎛⎜ ⎜ ⎜ ⎜⎝−M−1L00000MAAT0MAATS|A0|A000