{'fraction_non_alphanumeric': 0.04373518146999818, 'fraction_numerical': 0.013350355644720044, 'mean_word_length': 4.225080998665904, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Structure-Preserving Model Order Reduction for Index Two Port-Hamiltonian Descriptor Systems

We present a new optimization-based structure-preserving model order reduction (MOR) method for port-Hamiltonian descriptor systems (pH-DAEs) with differentiation index two. Our method is based on a novel parameterization that allows us to represent any linear time-invariant pH-DAE with a minimal number of parameters, which makes it well-suited to model reduction. We propose two algorithms which directly optimize the parameters of a reduced model to approximate a given large-scale model with respect to either the H-infinity or the H-2 norm. This approach has several benefits. Our parameterization ensures that the reduced model is again a pH-DAE system and enables a compact representation of the algebraic part of the large-scale model, which in projection-based methods often requires a more involved treatment. The direct optimization is entirely based on transfer function evaluations of the large-scale model and is therefore independent of the system matrices' structure. Numerical experiments are conducted to illustrate the high accuracy and small reduced model orders in comparison to other structure-preserving MOR methods.


Introduction
In this article, we generalize previously developed optimization-based structure-preserving model order reduction (MOR) algorithms for port-Hamiltonian differential algebraic systems (pH-DAEs) of differentiation index one [1] to pH-DAEs with index two, which occur very frequently in applications, see, e. g., [2]. To accommodate the hidden constraints that may occur in higher-index pH-DAEs, we adapt both the parameterization of the reducedorder model (ROM) and the optimization strategy compared to the case of index one pH-DAEs, as considered in [1].
As a system class, we consider linear constant coefficient pH-DAE systems of the form The research by P. Schwerdtner where E, J, R ∈ R n×n , G, P ∈ R n×m , S, N ∈ R m×m and the quadratic Hamiltonian is H(x) = 1 2 x T Ex. The following conditions are imposed on the system matrices: (i) The matrices J and N are skew-symmetric.
(ii) The passivity matrix W := R P P T S and the energy matrix E are symmetric positive semi-definite.
Note that we consider systems with Q = I n rather than the more general definition of pH-DAEs in [3] governed by a quadratic Hamiltonian of the form H(x) = 1 2 x T Q T Ex with Q T E ≥ 0. However, this does not pose any real restriction, since it has been shown that every linear timeinvariant pH-DAE may be reformulated to display this form; see [4] for a detailed discussion.
In this paper, we treat systems that satisfy the following additional assumptions: (i) the pencil sE − (J − R) is regular, i. e., det(sE − (J − R)) is not the zero polynomial; and (ii) all finite eigenvalues of sE−(J −R) have a negative real part.
PH-DAE systems of this type occur in energy-based modeling in many different physical domains, such as electrical systems [5], flow-problems [6,7], or mechanical multi-body systems [2,Example 12]. The pH formulation is particularly useful where systems from different physical domains or of different scales and accuracy are coupled. Since the pH structure is preserved for parallel and negative feedback interconnections of pH systems and the index of pH-DAEs is at most two as shown in [8], pH-DAEs can be safely coupled without increasing the index. In contrast, for general DAEs, the index may increase arbitrarily in such a coupling, and eigenvalues may move to the right-half plane, leading to an unstable system behavior. This emphasizes the usefulness of pH-DAEs when modeling subsystems in large networks of dynamical systems as in [6,9,10].
When a pH-DAE describes a complex physical process or a complex geometry arises during modeling, the system typically has a high state-space dimension n, which renders the direct simulation or model-based controller design computationally expensive. In this case, MOR is used to compute reduced-order models (ROMs) that lead to similar output trajectories as their full-order model (FOM) counterparts when the same inputs are provided. The ROM can then be used in simulations and for control instead of the FOM.
The input-to-output behavior of a linear system in the frequency domain is described by its transfer function. Under the assumption of regularity of the pencil sE −(J −R), the transfer function of the system in (1) is given by The typical use-case of pH systems in networks requires a structure-preserving MOR of pH-DAEs, i. e., the ROM must also be a pH system. Otherwise, the intuitive coupling to the rest of the network system is not possible at the ROM level. In particular, the ROM is supposed to be of the form where the matrices E r , J r , R r ∈ R r×r , G r , P r ∈ R r×m , S r , N r ∈ R m×m satisfy the same assumptions as the system in (1) with r n and for each admissible input function u we have that y ≈ y r . The latter condition can be expressed as the requirement that the transfer function of the ROM should be a good approximation of the transfer function of the FOM.
The transfer function H can be decomposed as where H sp is a strictly proper rational matrix function and H pol is a matrix polynomial of degree ν − 1. Here, ν denotes the differentiation index of the uncontrolled DAE. The additional polynomial part in a FOM transfer function poses an extra challenge for MOR of DAEs because it must be matched exactly in the ROM transfer function. Otherwise the error between FOM and ROM can become unbounded (see Section 2 for error measures). Therefore, in projection-based MOR, the algebraic part of the DAE is typically treated separately, see [11,12,13,14,15]. While for ν = 1, the additional polynomial part is always constant, such that the FOM system can still be approximated with an ODE system with feedthrough term, for ν = 2, the transfer function of a pH-DAE model includes a linear term, which requires a reduced-order DAE model. While the MOR of general higher-index DAEs has been extensively studied since the early 2000s (see [16] for balancing-based methods and [13,17] for interpolationbased methods), the structure-preserving MOR of pH-DAEs has still only partially been resolved. The methods in [11,18] require a separate treatment of the algebraic part in order not to destroy the pH structure and often lead to an accuracy several orders of magnitude lower compared to non-structure-preserving MOR methods. Furthermore, the method in [11] does not guarantee preservation of the pH structure when the FOM transfer function is improper and is only generally applicable to pH-DAEs where the hidden constraints are uncontrolled.
An alternative approach to pH structure-preserving MOR is passivity-preserving MOR, since the pH structure can often be recovered from a passive ROM [2]. However, passivity-preserving methods such as positive-real balanced truncation (PRBT), as in [19], often require an extra treatment of the algebraic equations which may impose an additional computational burden. Furthermore, transformations to the FOM matrices are needed that either require a priori knowledge of the given DAE structure or the construction of spectral projectors [15].
This article describes an optimization-based MOR approach for structure-preserving MOR that The paper is organized as follows. In Section 2 we provide a short background on MOR and explain the properties and different types of pH-DAEs using a staircase form. In Section 3 we present and explain our new method, which includes a parameterization of pH-DAEs and two optimization strategies implemented in Algorithms 1 and 2. Numerical experiments that illustrate the effectiveness of our approach are presented in Section 4.

Preliminaries
In this section we present some preliminary results for pH-DAEs that form the basis of the construction of the MOR algorithms. We then summarize several MOR techniques for DAEs and their properties.

Staircase form
Our first result resembles a staircase form for pH-DAEs, which was derived for general linear pH-DAE systems without input and output in [20] and for the system with input and output in [11], see also [3]. The staircase form allows us to determine the differentiation index of a given pH-DAE but is based on subsequent rank decisions which may be a challenging task in finite precision arithmetic. Fortunately, the staircase form often emerges naturally from the modeling process [2,21] or can be constructed using only a few permutations of the given system matrices. We present the real-valued version of the form as derived in [20].

Lemma 1. [20]
Consider a regular pH-DAE system as in (1). Then there exists an orthogonal matrix T ∈ R n×n such that is positive definite and the matrices J 41 and J 33 − R 33 are invertible (if the blocks are nonempty).
The transformed system is again a pH-DAE system with accordingly partitioned state vector

Model order reduction of higher-index systems
MOR for descriptor systems of possibly higher index have been proposed by many authors, see, e. g., [13,15,16,22,23,24,25], to mention but a few. In such methods, the error between the transfer function of the FOM and the ROM are typically measured in norms defined in the Hardy spaces RH p×m ∞ and RH p×m 2 of all proper (resp. strictly proper) real-rational p × m matrices without poles in the closed complex half-plane C + := {s ∈ C | Re(s) ≥ 0}. The typical norms are respectively, see, e. g., [26] for details. In many MOR methods one can bound or even optimize the transfer function error in these norms.
The iterative rational Krylov algorithm (IRKA) aims at minimizing the H 2 error, see [27]. A structure-preserving variant known as IRKA-PH, was derived in [28], which is extended to DAEs in [11] using the staircase form in Lemma 1. The strictly proper part and the Markov parameters of the original transfer functions are identified, such that the standard version of IRKA-PH can be applied. These algorithms are based on the tangential interpolation of the transfer function at iteratively defined interpolation points and tangential directions and the construction of projection spaces using an orthonormalization procedure based on rational Krylov subspaces, see [29] for an exhaustive discussion of such interpolation methods. Passivity-preserving model reduction methods tailored to DAEs were proposed in [11,19]. Positive-real balanced truncation (PRBT) [19] is based on the computation of the minimal solutions of two algebraic Riccati equations (or Lur'e equations). These minimal solutions take the role of the Gramians in classical balanced truncation (BT) model reduction, and the balancing and truncation procedure can be performed as in BT. While this method also admits an a priori gap metric error bound [30], it generally does not preserve the pH form.

Our method
Our method is based on optimizing the parameters (i. e., the matrix elements) of a low-order pH-DAE such that its transfer function matches the transfer function of a given large-scale pH-DAE a best as possible. As error measures we consider both the H 2 and the H ∞ norms. For this we need to (i) determine a parameterization that encompasses all the pH-DAE types that are described in Lemma 1, (ii) impose a matching of the polynomial part of the given transfer function to keep the errors well-defined and bounded, and (iii) define update rules for the parameters such that H 2 or H ∞ errors are minimized.
We set up a parameterization that allows the strictly proper part of the transfer function to be tuned independently of the polynomial part in Section 3.1. This way, matching of the polynomial part and minimizing the errors can be addressed separately. We present several alternative methods of computing the polynomial part of a potentially large-scale pH-DAE in Section 3.2. Methods of both H ∞ and H 2 approximation are then given in Sections 3.3 and 3.4, respectively.

Parameterization
Before we set up the parameterization, we shall first recall a few properties of transfer functions of passive systems which are positive real. A real-rational matrix H is referred to as positive real, if (iii) For each ω ∈ R for which iω is not a pole of H, it holds that Lemma 2. Let a real-rational positive real m × m matrix H be given. We can decompose it into the sum Proof. We refer to [32] for the proof of (a). To show (b), we decompose H into its proper and improper parts, and for each ω ∈ R for which iω is not a pole of H we obtain where we have used (a). Consequently, the proper part H p is positive real.
Remark 1. If all finite eigenvalues of the pencil sE − (J − R) have negative real part as imposed by our assumptions, then its transfer function H cannot have finite poles on the imaginary axis, in particular, we have This result allows us to consider the improper part D 1 · s separately from the proper part of the given transfer function, as the proper part is still positive real and can thus be parameterized separately. For pH-DAE systems as in (1), we can even make further assertions about the structure of the proper part H p . Theorem 1. The proper part H p of the transfer function of any pH-DAE can be realized by an implicit ODE system of the form with x p : R → R n2 and E p > 0, where n 2 is the corresponding dimension as in Lemma 1.
Proof. The claim that the proper part of any transfer function of a pH-DAE admits a realization of the form (3) follows immediately, as any proper positive real real-rational matrix is realizable by a passive ODE system, which can in turn be transformed into pH form [33]. In the Appendix, we show that we can build a realization with state-space dimension n 2 by deriving it directly from the staircase form in Lemma 1.
In Theorem 2 we exploit the presented properties to derive a novel approach to parameterizing pH-DAEs of differentiation index two. Therein, the functions vtu(·) and vtsu(·) map vectors row-wise to appropriately sized upper triangular and strictly upper triangular matrices, respectively. The function names are abbreviations for vector-toupper and vector-to-strictly-upper, respectively. The function vtf m (·) is a simple reshaping operation, which maps a vector in R n·m to a matrix in R n×m . A detailed description of these functions can be found in [34, Definition 3.1].
Theorem 2. Let θ ∈ R n θ be a parameter vector parti- Then the parameter-dependent DAE satisfies the pH structure conditions. Its transfer function H r is given by where H p,r (s, θ) denotes the transfer function of the proper pH-ODE system Σ p,r (θ) : Proof. For ease of notation, we omit the argument θ in the system matrices of (4). The fact that (5) is a pH-DAE system of index two follows directly from the composition of the reduced-order matrices and the parameterization in (4). The system may be reformulated aṡ y p,r (t) = (G p,r + P p,r ) T x 1 (t) where we have used the fact that x 3 (t) = −ẋ 2 (t) = L Tu (t) (assuming that u is differentiable). This leads to the transfer function H r (s) = H p,r (s) + LL T · s, with the proper part H p,r (s) = (G p,r + P p,r ) T (sI r − (J p,r − R p,r )) −1 (G p,r − P p,r ) The assertion follows from the fact that H p,r is the transfer function of the system in (6), which clearly fulfills the pH structure conditions.

Remark 2.
We highlight that this parameterization naturally carries over to proper pH-DAEs, i. e., with either index ν < 2 or G 4 = 0. For these systems we have that = rank(L) = 0 and therefore it is sufficient to parameterize the system with an implicit pH-ODE as in (6). This case is discussed in more detail in [1].
Remark 3. Note that for the proper subsystem of the reduced model Σ p,r (θ), we assume that E p,r (θ) = I r . This is not a restriction, since every pH-ODE system may be transformed into such a form using, for instance, the Cholesky factor of E p,r (θ).

Computation of D 0 and D 1
To keep the H ∞ error between the FOM and ROM transfer functions bounded, the polynomial coefficient D 1 of the FOM transfer function must be matched exactly in the ROM. To obtain a bounded H 2 error, both D 0 and D 1 must be matched exactly. Since our parameterization allows the polynomial coefficients to be assigned and the remaining free parameters to be independently optimized, it remains to compute D 0 and D 1 efficiently for the FOM transfer function.
One approach is to determine D 0 and D 1 by the method outlined in [15], where these coefficients can be read off a block-Hankel matrix constructed using the solutions of two projected discrete-time Lyapunov equations. Another approach consists in transforming the FOM to the almost Kronecker canonical form derived in [20], see also Lemma 3 in the Appendix.
Alternative approaches determine D 0 , D 1 by sampling the original transfer function H; hence these are independent of the specific representation of the FOM. Assume that two distinct imaginary sampling points s 1 = iω 1 , s 2 = iω 2 with two sufficiently large values ω 1 = ω 2 are given. As shown in [35], we have that which follows from lim s→∞ H sp (s) = 0. We can then obtain an estimate for D 1 as Similarly, we have that which yields an estimate for D 0 given by The estimation quality generally depends on the choice of ω 1 , ω 2 , and we can adapt these sampling points in an iterative manner as in [36], until a certain tolerance is met. Another way of enhancing the accuracy is to include more than two sampling points, using the Loewner framework [35].

H ∞ approximation
Since we assume that the pencil sE − (J − R) is regular and its finite eigenvalues are in the open left half plane, we can proceed as in [34] and obtain a good H ∞ approximation by minimizing for decreasing values of γ > 0 the function L(θ; H, γ, S) with respect to θ, where Here, S ⊂ iR is a set of sample points at which the original and reduced transfer functions are evaluated and σ j (·) denotes the j-th singular value of its matrix argument. The procedure for the H ∞ approximation pH-DAEs of index two is described in Algorithm 1. It is based on repetitively minimizing L in conjunction with a bisection over γ. Note that L(·; H, γ, S) attains its global minimum at zero when the error H(s i ) − H r (s i ) 2 at all sample points s i ∈ S is smaller than γ. With appropriately chosen sample points (an adaptive sampling procedure is proposed in [37]), this in turn is an indication of H − H r (·, θ) H∞ < γ. Thus, finding the minimal γ such that L(·; H, γ, S) can be minimized to zero (as in Algorithm 1) is a reasonable strategy for determining a ROM with a small H ∞ error. In [34], we discuss the further benefits of this general approach, also in comparison to a direct minimization of the H ∞ norm.

H 2 approximation
For H 2 -optimal model reduction, we first have to ensure that H − H r (·, θ) ∈ RH m×m 2 (see Section 3.2). In addition to the condition that L(θ)L(θ) T = D 1 , we therefore also require that This may be achieved by fixing θ L , θ N and m(m+1)/2 parameters in θ W once D 0 and D 1 have been obtained using any of the methods in Section 3.2 (see [1] for details). The remaining n θ − m(m + ) parameters in θ can be tuned to minimize the error H − H r (·, θ) H2 , as initially proposed in [38]. Assume that the eigenvalues of J p,r (θ)−R p,r (θ) are simple, and consider the spectral decomposition (J p,r (θ) − R p,r (θ))Z(θ) = Z(θ)Λ(θ), where Λ(θ) = diag(λ 1 (θ), . . . , λ r (θ)) and Z(θ) contains the right eigenvectors as columns. The transfer function H r (·, θ) may then be represented by the partial fraction and where e i denotes the i-th standard basis vector of R r . As shown in [39, Theorem 2.1], we may then express the H 2 error by where H sp again denotes the strictly proper part of H.
Since H sp 2 H2 does not depend on θ, it can be neglected in the optimization. Consequently, we define the objective functional which can be evaluated efficiently, since it only requires the solution of the reduced-order eigenvalue problem in (8) and r evaluations of H sp at −λ i (θ), i = 1, . . . , r. Note that to evaluate H sp , we do not require a state-space representation of the strictly proper part. Instead, we indirectly evaluate H sp by for all s ∈ C. Since the original matrices are sparse, this is much more efficient than an explicit computation of the strictly proper part of the FOM, as shown, for example, in [20]. The objective functional F(·; H) can then be optimized with respect to θ using gradient-based approaches. For the derivation of the gradient and possible initialization strategies, we refer to [1] and summarize the procedure in Algorithm 2.

Numerical examples
In this section, we illustrate the properties of our new methods using several numerical tests based on well-known benchmarks.

Benchmark models
The first test case we consider is the instationary flow of incompressible fluids on the spatial domain Ω = (0, 1) 2 Algorithm 2: PROPT-H 2 Input : FOM transfer function H, reduced order r ∈ N. Output: Reduced-order model as in Theorem 2 with parameter θ fin .
4 Construct the ROM with θ fin as in Theorem 2.
with velocity vector v, pressure p, viscosity ν > 0, external forces f and a convective term with driving velocity a. As shown in [18], spatial semi-discretization of the Oseen equations by a finite difference method with associated non-slip boundary conditions and initial velocity v 0 ∈ R 2 , i. e., leads to a pH-DAE of index two and with n 3 = 0 in the staircase form of Lemma 1. Additionally, we have that G 4 = 0, which makes the transfer function H strictly proper. We consider two models: Oseen-1 with n = 279 and Oseen-2 with n = 7399. PH-DAEs of index two also occur in the network modeling of RCL circuits by directed graphs, as described in [40]. The second test case we consider relates to RCL ladder networks as shown in Figure 1, which are part of the software package PortHamiltonianBenchmarkSystems 1 . Here,n denotes the number of loops in the system. If we choose  Figure 3 shows the H ∞ and H 2 errors for the Oseen-1 model and both RCL systems. To compare our methods with the IRKA-PH and PRBT algorithms, we plot these errors over the dimension r of the proper part of the ROM. We highlight that the ROMs of our methods are always minimal realizations with state-space dimension r +2 , see Theorem 2. In the case of the PRBT algorithm, minimal realizations can also be computed but may require the solution of discrete-time projected Lyapunov equations, see [19]. Also note that the ROMs produced by PRBT are generally not in pH form.

Results
Let us first consider the MOR results on the strictly proper Oseen-1 model, shown in Figures 3 (a) and (b). It can be observed that both PROPT-H 2 and SOBMOR-H ∞ lead to more accurate models in both the H 2 and the H ∞ norms compared to IRKA-PH, especially when the reduced model order is increasing. Compared with PRBT, SOBMOR-H ∞ leads to slightly smaller H ∞ errors except for r = 10, while PROPT-H 2 leads to slightly smaller H 2 errors up to reduced orders of 8. Figures 3 (c) and (d) show the results for the RCL models. While the difference between the improper FOM and improper ROM transfer function is itself proper, since we ensure a matching of the polynomial parts, the computation of the H 2 norm using the Matlab function h2norm fails. Therefore, we only report the H ∞ errors. In terms of the H ∞ norm, the ROMs computed with SOBMOR-H ∞ have the smallest errors for both RCL models and all reduced model orders. The H ∞ error of the PROPT-H 2 ROMs stagnates after a reduced model order of 16 and is worse in general. However, this is expected, since PROPT-H 2 aims at minimizing the H 2 error. Figure 4 shows the results for the Oseen-2 model. Since no implementation of PRBT that exploits sparsity is currently publicly available, we only report results obtained with our methods and IRKA-PH. Our findings are comparable to the smaller Oseen-1 model: The error transfer functions of IRKA-PH and PROPT-H 2 are similar for the small reduced model order r = 5, while for r = 10, the PROPT-H 2 error is clearly smaller than the error of the IRKA-PH model across a wide frequency range. SOBMOR-H ∞ aims at minimizing the maximum error across all frequencies and it can be observed that the maximum errors of the SOBMOR-H ∞ ROMs are well below the maximum errors of the ROMs obtained with the other methods for both reduced orders.

Conclusion
We have presented a flexible MOR approach for H 2 and H ∞ approximation of higher index pH-DAEs. Our approach is based on a novel parametrization that can provide a pH realization for any pH descriptor system with an efficient representation of the algebraic part. An adaptation of previously created optimization-based MOR methods allows for the approximation of potentially improper transfer functions. A comparison to state-of-the-art methods shows that our optimization-based approach leads to highly accurate ROMs that are minimal realizations and guaranteed to fulfill the pH structural constraints.
The matrices are partitioned in the same way as in Lemma 1 and, if present, the diagonal block matriceš E 11 ,Ě 22 are symmetric positive definite.
Proof of Theorem 1. Considering the transformed system (10) and a block diagonalization of sĚ −Ǎ yield the transfer function which reveals the split into the proper and improper parts, respectively. From the definition of T 1 , T 2 in [20] we have thatĚ which proves that D 1 =Č 4Ě11B4 = D T 1 ≥ 0. For n 2 = 0, the claim follows immediately from H(0) + H(0) H = D 0 + D T 0 ≥ 0. Now let us first consider the most general case where n 2 , n 3 > 0. The remaining block matrices from (10) are then given by which again follows from the definition of T 1 , T 2 in [20]. Define the matrices Γ, W ∈ R (n+m)×(n+m) and let the matrices be partitioned as in Lemma 1. Our proof is based on the observation that This is a natural generalization of a similar observation for linear dissipative Hamiltonian systems (see [21]) to port-Hamiltonian systems with power-collocated input-output pairs. We will now show that such a decomposition into skew-symmetric and symmetric positive semidefinite parts not only exists for the proper subsystem as well but may be obtained by structure-preserving manipulations of the sum Γ + W . At first, let P π ∈ R (n+m)×(n+m) define a permutation matrix which permutes the third and fifth block rows and columns in the sum Γ + W such that Ψ = P T π Γ + W P π = Ψ uu Ψ ul Ψ lu Ψ ll , with Ψ ll = − J 33 + R 33 = − A 33 . Since the matrix A 33 is nonsingular, we may block-diagonalize Ψ with invertible matrices Note that Ξ still has a positive semidefinite symmetric part since the Schur complement preserves this property [21,Corollary 4.3]. Finally, it is easy to show that we may compute the proper system matrices via a transformation of Ξ with the full-rank matrix U ∈ R (n+m)×(n2+m) such that Hence, we obtain the proper system matrices by a series of permutations, block-diagonalization via Schur complements and congruence transformations of Γ + W . Since each of these manipulations preserves the positive semidefiniteness of the symmetric part, we obtain a port-Hamiltonian representation of the proper subsystem via The fact that E p =Ě 22 > 0 proves the claim for n 3 > 0. Similar arguments apply to pH-DAEs of differentiation index two for which n 3 = 0 and where we have that This concludes the proof.