Model Reduction of Linear Differential–Algebraic Equation Systems Using Constructed Spectral Projectors

This contribution introduces linear differential–algebraic equation (DAE) systems and provides the explicit construction steps of spectral projectors for DAEs with indexes 1 or 2. Furthermore, an interpolatory projection-based model order reduction method using spectral projectors is applied to an eddy current problem as well as an electrical circuit to compute a reduced order model.


I. INTRODUCTION
T HE input-output behavior of many dynamic linear models can be described by using linear differential-algebraic equation (DAE) systems, for example, semi-discretized partial differential equations (PDEs) or electrical circuits using a modified nodal analysis (MNA). Here, in the case of passive systems, the algebraic part results in a dominating linear behavior of the transfer function in the upper frequency range. However, some model order reduction (MOR) techniques are unsuitable to approximate the behavior of the algebraic part of the transfer function, e.g., the rational Krylov method [5]. Therefore, different MOR techniques for solving DAE systems exist that require spectral projectors to decompose the DAE into a differential and an algebraic part, e.g., [2], [5]. However, Gugercin et al. [5] use solely analytical constructed spectral projectors, which restrict the introduced MOR method to DAE systems with special structures. In addition, the presented approach in [5] preserves the algebraic part, which results only in a partial degrees of freedom (DOFs) reduction.
The novelty of this contribution consists in the usage of spectral projectors constructed according to [1] and the application of the H 2 optimal interpolatory MOR [5] to the differential part of the DAE system. Furthermore, the algebraic part is reduced to increase the compression level.
In the following, Section II introduces DAE systems for single-input single-output (SISO) as well as the definition of DAE index and spectral projectors. In the case of DAE systems with index 1 or 2, Section III demonstrates the explicit construction of spectral projectors based on [1]. Furthermore, Section IV introduces the H 2 optimal interpolatory MOR technique using spectral projectors. This proposed MOR method is applied to an index 1 and 2 DAE system to compute a reduced order model (ROM), where the results are presented in Section V. Finally, Section VI provides a concluding discussion of the proposed approach.
with matrices E, A ∈ R N ×N , vectors b, c ∈ R N , d ∈ R, as well as time-dependent input u(t) ∈ R, solution x(t) ∈ R N , and output y(t) ∈ R, respectively. If there exists a λ 0 ∈ C, such that det(λ 0 E−A) ̸ = 0, then the matrix pencil (E, A) is denoted as regular. In this case, there exist nonsingular matrices S, T ∈ R N ×N , such that E, A can be represented by the Weierstraß-Kronecker canonical form [5] where n f + n ∞ = N , I n f ∈ R n f ×n f , and I n ∞ ∈ R n ∞ ×n ∞ are identity matrices, and J ∈ R n f ×n f , N ∈ R n ∞ ×n ∞ are matrices in Jordan canonical form. Moreover, N is a nilpotent matrix with index µ ≤ n ∞ , i.e., N, . . . , N µ−1 ̸ = 0 and N µ = 0. Here, µ is called the (Kronecker) index of (1) whereby, e.g., many semidiscretized parabolic PDEs and electrical circuits using MNA result in µ ≤ 2. Further, n f and n ∞ denotes the number of generalized finite and infinite eigenvalues of the regular matrix pencil (E, A), respectively. Using Laplace transform with s ∈ C, the transfer function where the strictly proper (sp) part H sp (s) corresponds to an ordinary differential equation (ODE), and the improper (imp) part H imp (s) is a polynomial of degree at most µ [1]. This decomposition can be implemented effectively using left (ℓ) and right (r ) spectral projectors III. CONSTRUCTION OF SPECTRAL PROJECTORS As described in [1], spectral projectors P ℓ , P r can be iteratively constructed using canonical projectors. The construction depends on the index µ of (1), and it is summarized in the following for µ ∈ {1, 2}. However, [1] provides the general construction procedure.
This construction leads to singular matrices E 0 , . . . , E µ−1 and a nonsingular matrix E µ according to [1,Theorem 1]. Therefore, [2] provides an efficient approach to compute the projector Q k onto ker E k , which is summarized in Algorithm 1 using MATLAB 1 syntax and functions. The essential step in the computation of P r is to construct canonical projectors Q 0,c and Q 1,c onto ker E 0 and ker E 1 , respectively, fulfilling (4) and other properties detailed in [1]. 1) Case µ = 1 [1]: For the canonical projector Q 0,c = Q 0 E −1 1 A 0 , the right spectral projector is given due to P r = (I N − Q 0,c ).

B. Construction of Left Spectral Projector P ℓ
The construction of P ℓ for µ ∈ {1, 2} follows in the same way as Section III-A, but with the difference that the initial matrices of (4) are set to E 0 = E ⊤ and A 0 = A ⊤ .   (1), provides, for n ≪ N , the ROM where E n = W ⊤ EV, A n = W ⊤ AV, b n = W ⊤ b, and c ⊤ n = c ⊤ V using appropriate projection matrices V, W ∈ R N ×n . Furthermore, H n (s) = c ⊤ n (s E n − A n ) −1 b n + d defines the reduced transfer function. As a specific method, the H 2 optimal interpolatory MOR method [5] is chosen, which relates to a non-convex optimization problem using H 2 -norm, to find V, W for stable reduced order SISO models, such that where ∥H ∥ H 2 := ((1/2π ) ∞ −∞ |H ( jω)| 2 dω) (1/2) . Here, the used approach of the implemented H 2 optimal interpolatory MOR technique based on [5, Algorithm 4.1] and Algorithm 2 summarizes the main steps to compute the ROM using MAT-LAB syntax, functions and toolboxes sssMOR [6] and MOR-LAB [7]. Therefore, in accordance with Section III, spectral projectors P ℓ , P r are constructed in MATLAB to decompose the transfer function H (s) into H sp (s) and H imp (s). Since the strictly proper part can be associated with an ODE, the iterative rational Krylov algorithm (IRKA) is applicable to the strictly proper part to compute an ROM of H sp (s). Furthermore, the polynomial part of H (s) corresponding to H imp (s) can be approximated by an ROM using techniques for the generalized discrete-time dual Lyapunov equations [8]. Finally, according to Algorithm 2, matrices V, W ∈ R N ×n are provided to compute the ROM, where n = n sp + n imp is decomposed into the DOF of the strictly proper and improper part, respectively. The benefit of Algorithm 2 is the simplicity of IRKA, and if IRKA converges, then in numerous applications, a rapidly convergence has been observed [5]. But, IRKA has also some drawbacks, for example, it is not guaranteed to converge to a Fig. 1. Investigation of inductive coupling of two wire coils using 2-D finiteelement approach, which result in an index 1 DAE system. local minimum, and the H 2 norm does not decrease monotonically during the iteration. Furthermore, unstable ROMs may occur, since, according to [5], IRKA satisfies only the first-order necessary optimality condition of (6). In addition, the choice of initial values has a strong influence on the convergence behavior.

V. NUMERICAL RESULTS
In the following, the presented approach is applied to two test problems corresponding to an index 1 and 2 DAE system, respectively. Here, spectral projectors P ℓ , P r for µ ∈ {1, 2} are constructed according to Section III, and the H 2 optimal interpolatory MOR is applied in accordance with Section IV. Therefore, computational costs and runtimes are provided, whereby the computation was performed on a Windows machine equipped with an intel processor (i7-6500U CPU@2.50 GHz) and 16 GB RAM. Fig. 1 depicts the test problem of two parallel wire coils with a distance of b = 10 mm apart. Each wire has a diameter of D = 2 mm, and the outer radius of the wire coils is R = 30 mm. In the first wire coil, an impressed current density J i is prescribed corresponding to the single input u(t) = i 0 (t). The second short-circuited wire coil is implemented with a conductivity γ = 56 · 10 6 ( m) −1 . Here, the single output y(t) corresponds to the eddy current i 1 (t) = 1 γ (−∂ t A) • dS, which result in the transfer function H ( j 2π f ) = ((I 1 ( j 2π f ))/(I 0 ( j 2π f ))), where I 0 (s) and I 1 (s) correspond to the Laplace transform of i 0 (t) and i 1 (t). As noted in Fig. 1, this application is modeled with an index 1 DAE system using the 2-D magneto-quasistatic A * formulation taking advantage of the radial symmetric geometry, where the spatial finite-element semi-discretization is implemented with the first-order nodal elements using openCFS [9]. Moreover, the 16 initial expansion points s 0 ∈ {± j 2π 10 k : k ∈ N, 1 ≤ k ≤ 8} guarantee the construction of real-valued matrices V f , W f , V imp , and W imp . Fig. 2 depicts the strictly proper parts |H sp | and |H sp,n sp | and the reduced improper part |H imp,n imp |. Here, the full order improper part relates to n ∞ = 7814 DOF corresponding to the air domain in the test problem, which confirms the importance of an ROM for the improper part, since a reduction to n imp = 1 is possible. Fig. 3 shows that the ROM leads to an accurate  1 ( j 2π f ))/(I 0 ( j 2π f ))).

Fig. 3.
Transfer functions H ( j 2π f ) = ((I 1 ( j 2π f ))/(I 0 ( j 2π f ))), H n ( j 2π f ) as well as error relating to FOM and ROM. approximation of the FOM and reduces the DOF from N = 9440 to n = 17 significantly. Therefore, the construction of the spectral projectors P ℓ , P r needed a runtime of 189.12 s, and the total runtime of the ROM is 687.90 s, including 3.62 s for the strictly proper part. Furthermore, the computation of one frequency step of the transfer function requires a runtime of 224.61 ms for FOM and 2.01 ms for ROM using the sssMOR function "freqresp" [6].

B. Index 2 DAE System
Following [5], the electrical circuit according to Fig. 4 is used to investigate the H 2 optimal interpolatory MOR method regarding Section IV, but with the difference, that instead of analytical spectral projectors, P ℓ and P r are constructed numerically, and an ROM of the improper part is computed.
|H sp | and |H sp,n sp | using the 16 initial expansion points s 0 ∈ {± j 2π 10 k : k ∈ N, 1 ≤ k ≤ 8}, as well as the reduced improper part |H imp,n imp |. Here, the full order improper part relates to n ∞ = 502 DOF, which allows a reduction of the improper part to n imp = 2 DOF. Fig. 6 shows that the ROM leads to an accurate approximation of the FOM and reduces the DOF from N = 1502 to n = 18. Therefore, the construction of the spectral projectors P ℓ , P r needs a runtime of 3.35 s, and the total runtime of the ROM is 3.10 s, including 0.59 s for the strictly proper part. Furthermore, the computation of one frequency step of the transfer function requires a runtime of 10.14 ms for FOM and 1.88 ms for ROM using the sssMOR function "freqresp" [6].
VI. CONCLUSION Decomposing a DAE system into the algebraic and differential part leads to the possibility of constructing an ROM for both parts separately, which result in a more accurate ROM. Here, the availability of spectral projectors P ℓ and P r is essential for the practical implementation. Therefore, the presented approach for index 1 and 2 DAE systems, according to Section III, provides the usage of constructed spectral projectors P ℓ , P r with the major advantage, that it is independent of the DAE system structure. Furthermore, it provides the DAE index µ as an additional result, which is an important property, e.g., choosing an appropriate time stepping method solving a DAE. Therefore, according to Section III, Algorithm 1 can be used to determine µ by checking if rank E k = N .
On the basis of the results obtained, the discussed approach can be recommended for different applications resulting in a DAE system with less than 10 4 DOF. However, large-scale problems with more than 10 4 DOF require some improvements with respect to computational costs, e.g., using parallel programming or avoiding the computation of inverse matrices, whenever possible.
In future work, the proposed approach can be investigated for multi-input multi-output (MIMO) DAE systems, since [5, Algorithm 4.1] is not restricted to SISO DAE systems.