Hankel-Norm Approximation of Large-Scale Descriptor Systems

The Hankel-norm approximation is a model reduction method which provides the best approximation in the Hankel semi-norm. In this paper the computation of the optimal Hankel-norm approximation is generalized to the case of linear time-invariant continuous-time descriptor systems. An efficient algorithm is developed by refining the generalized balanced truncation square root method. For a wide practical usage, adaptations of the introduced algorithm towards stable computations and sparse systems are made as well as an approach for a projection-free algorithm. To show the approximation behavior of the introduced method, numerical examples are presented.


Introduction
Many different real-world applications, like chemical processes, electrical circuits and networks, or computational fluid dynamics, naturally lead to models, described by systems of differentialalgebraic equations. Since experiments can be very costly, time-consuming, and expensive, these models are used for simulations and the design of controllers. The modeling process often results in linear time-invariant continuous-time descriptor systems of the form with E, A ∈ R n×n , B ∈ R n×m , C ∈ R p×n , D ∈ R p×m . Here, u(t) ∈ R m are the inputs of the system, which influence the generalized states x(t) ∈ R n to get the desired outputs y(t) ∈ R p . Throughout this paper, it is assumed that the matrix pencil λE − A is regular, i.e., there exists at least one λ ∈ C such that det(λE − A) = 0. In this case, and with the initial condition Ex(0) = 0, the input-output behavior of the system (1) in the frequency domain can be described via the system's transfer function The quintuple (E, A, B, C, D), consisting of matrices from (1), defines a realization of (1) and its transfer function (2). Usually, the numbers of inputs and outputs are very small in contrast to the number of differential-algebraic equations and generalized states n, which quickly enlarges due to different reasons, e.g., the model shall provide a required accuracy. Because of that, the usage of complete models often reaches the limits of computational resources like memory and computation time. Since the acquired data for the model usually contain a huge amount of redundancies, it is possible to approximate the original model by a new system with a much smaller order. The task of model reduction is to construct a reduced-order descriptor system Eẋ(t) =Âx(t) +Bu(t), of order r ≪ n, such that the input-output behavior of the original system (1) is approximated. Many model reduction techniques were originally developed for the standard system case, where the descriptor term E is the identity matrix I n (or at least nonsingular). But in recent years, quite a few of these methods have been extended to the case of descriptor systems with singular E matrices. There are different approaches for the construction of (3), e.g., matrix equations can be used to determine a measure for truncatable states [6], or the transfer function can be approximated by rational interpolation [12]. A special technique of model reduction is the computation of the optimal Hankel-norm approximation (HNA). This technique actually provides a best approximation in the Hankel semi-norm. Based on the work of Adamjan, Arov, and Krein about the approximation of Hankel matrices [1], an algorithm for the computation of the HNA for standard systems was introduced by Glover in [11].
A generalization of the HNA to the descriptor system case was already mentioned by Cao, Saltik, and Weiland in [10]. They are using the Weierstrass canonical form for an explicit construction of reduced decoupled subsystems. The main problem of this method is the computation of the Weierstrass canonical form which is numerically costly and unstable. Also additional conditions, like C-controllability and C-observability of the system, have to be assumed.
In this paper, a new efficient algorithm for the computation of the generalized Hankel-norm approximation (GHNA) will be proposed. Our main contributions are twofold: 1. We generalize the concept of all-pass transfer functions to descriptor systems (Theorem 1).
2. We derive new and reliable numerical implementations of the GHNA that also allow the application of the Hankel-norm approximation method to large-scale problems with sparse coefficient matrices as they arise, e.g., from systems with dynamics described by semi-discretized unsteady partial differential equations.
Therefor, in Section 2 the mathematical background of linear descriptor systems is recalled. Then, the HNA method for the standard system case is introduced in the first part of Section 3. Afterwards, the generalized balanced truncation is reviewed and used for the construction of the new GHNA method. The numerical difficulties and adjustments are discussed in Section 4 for usable implementations of the method. Two different implementations of the method are then tested on numerical examples in Section 5. In Section 6, the conclusions of this paper can be found.

Mathematical Basics
For regular matrix pencils λE−A, the Weierstrass canonical form always exists: there are invertible matrices W, T ∈ C n×n such that where J and N are both in Jordan canonical form, J is regular, and N is nilpotent with index ν.
The numbers n f and n ∞ are the dimensions of the deflating subspaces corresponding to the finite and infinite eigenvalues of λE − A, respectively. Then, the spectral projectors onto the left and right deflating subspaces corresponding to the finite eigenvalues of the matrix pencil λE − A can be defined as with W and T from the Weierstrass canonical form (4). Another necessary assumption is the c-stability of the matrix pencil λE − A, i.e., the matrix pencil λE − A is regular and all finite eigenvalues of λE − A lie in the open left half-plane. In this case, the proper controllability and observability Gramians are defined as the unique, positive semidefinite solutions of the projected generalized continuous-time Lyapunov equations with P ℓ and P r the spectral projectors corresponding to the finite eigenvalues (5); see [18]. Furthermore, the improper controllability and observability Gramians are given as the unique, positive semidefinite solutions of the projected generalized discrete-time Lyapunov equations with Q ℓ = I n − P ℓ and Q r = I n − P r the spectral projectors onto the left and right deflating subspaces corresponding to the infinite eigenvalues of the matrix pencil λE − A; see [18].
Using the system Gramians, the set of Hankel singular values is defined in the following; see [14].
Definition 1. The square roots of the n f largest eigenvalues of G pc E T G po E denoted by ς 1 ≥ ς 2 ≥ · · · ≥ ς n f are the proper Hankel singular values of (1). The square roots of the n ∞ largest eigenvalues of G ic A T G io A denoted by θ 1 ≥ θ 2 ≥ · · · ≥ θ n∞ are the improper Hankel singular values of (1).
In case of a non-singular descriptor term E, the proper Hankel singular values are the classical Hankel singular values of the system. Therefor, an equivalent energy interpretation of the proper Hankel singular values exists which proposes the truncation of states corresponding to small proper Hankel singular values, which are difficult to control and observe. Unfortunately, this does not hold for the improper Hankel singular values since these correspond to the constraints of the system. The truncation of non-zero improper Hankel singular values may result in physically meaningless systems.
There exist diverse concepts of controllability and observability for descriptor systems. For this paper, we restrict ourselves to the following definitions; see, e.g., [18]. 2. C-controllable if the system is R-controllable and rank E, B = n.
4. C-observable if the system is R-observable and The relation between these controllability, observability notions and the system Gramians is given in [19,Theorem 2.3]. Especially, all proper Hankel singular values are non-zero if and only if the system is R-controllable and R-observable.
The mapping from past inputs u − : (−∞, 0] → R m to future outputs y + : (0, +∞] → R p is described by the Hankel operator y + = Hu − . A generalization of this operator to the case of descriptor systems can be found in [10]. The measure of the influence of past inputs on future outputs in the L 2 -norm leads to the definition of the Hankel semi-norm for descriptor systems. Definition 3. The Hankel semi-norm of a system G is given by where W ν−1 2 (−∞, 0] denotes the Sobolev space of ν − 1 times weakly differentiable functions w.r.t. the L 2 inner product on the interval (−∞, 0] and . L2 the L 2 -norm. In case of an invertible descriptor term E, the Hankel semi-norm (10) simplifies to where ς max (G) is the largest Hankel singular value of the system G.

Algorithm for Standard Systems
First, the algorithm for the standard system case, introduced by Glover in [11], is considered. Therefor, a balanced minimal realization of the given standard system (I nmin , A, B, C, D) is assumed, where n min is the McMillan degree of the system, i.e., the order of its minimal realization. The computation is usually done by the balanced truncation square root method. Since the resulting system is balanced and minimal, the system Gramians are equal and diagonal G pc = G po = diag(ς 1 , ς 2 , . . . , ς nmin ), with ς 1 , . . . , ς nmin all non-zero Hankel singular values of the system. Next, the system is partitioned by the order r such that with k ≥ 1 being the multiplicity of the (r + 1)-st Hankel singular value. The Gramians are reordered to separate the block with the (r + 1)-st Hankel singular value aš withΣ = diag(ς 1 , . . . , ς r , ς r+k+2 , . . . , ς nmin ). Accordingly to (11), the remaining system matrices have to be permuted and partitioneď Then, the partitioned system is transformed by the following formulas with Γ =Σ 2 −ς 2 r+1 I nmin−k and U = (C T 2 ) † B 2 . Here, M † denotes the Moore-Penrose pseudo-inverse of a matrix M . This system is constructed such that the error transfer function E = G −G is scaled all-pass withG the transfer function of (12), i.e., it holds for all s ∈ C that are not poles of E(s) or E T (−s). In this case, the approximation error satisfies The transfer functionG of (12) has exactly n min − k − r unstable poles. As last step, an additive decomposition ofG is computed such thatG = G h + G + , where G + is the anti-stable part of order n min − k − r and G h is the stable part of order r. Since the Hankel semi-norm only depends on the stable part of the system, the error (14) in the Hankel semi-norm does not change if the unstable part is removed, such that

Computing a Balanced Realization for Descriptor Systems
As for the standard system case, for descriptor systems, a balanced conditionally minimal realization is needed. The term "conditionally" minimal means that the order of the system is minimal except of the reduction of the index-1 parts in E, see [17]. The computation is done using the generalized balanced truncation square root method (GBT(SR)). The basic idea of this method is the computation of a balanced realization and the truncation of unnecessary states.  [14]. The number of non-zero improper Hankel singular values is equal to the rank of the matrix G ic A T G io A, which can in fact be bounded by with ν, index of the system, m, number of inputs, p, number of outputs, and n ∞ , dimension of the deflating subspace corresponding to the infinite eigenvalues of λE − A. So for large n ∞ and usually small ν, the descriptor system (1) can be reduced significantly by the truncation of zero improper Hankel singular values. One method to compute the balanced truncation of a descriptor system is the square root method. Therefor, consider the skinny singular value decompositions 3 have orthonormal columns and the diagonal matrices Σ 1 , Σ 2 , and Θ 3 contain the non-zero proper and improper Hankel singular values, respectively. The partition of the proper Hankel singular values is chosen such that Σ 1 contains all the desired Hankel singular values and Σ 2 the undesired ones. By using the singular value decompositions in (17) and (18), the following transformation matrices can be defined where ℓ = ℓ f + ℓ ∞ is the sum of the number of desired proper Hankel singular values ℓ f and the non-zero improper Hankel singular values ℓ ∞ . The transformed realization is of order ℓ and balanced with the set of Hankel singular values contained in Σ 1 and Θ 3 . The resulting matrix pencil λÊ −Â is a resembling of the Weierstrass canonical form (4), such that Due to the reason that only the zero improper Hankel singular values have been truncated, the polynomial part of the system G has not changed. So it can be shown that the same error bound as for the classical balanced truncation method holds. LetĜ be the transformed descriptor system (20), then it holds with ς k (G) the k-th proper Hankel singular value of G.

Hankel-Norm Approximation of Descriptor Systems
As for the standard case, the GHNA method for descriptor systems is based on the construction of an error system with all-pass transfer function (13). The following theorem provides an algebraic characterization of descriptor systems with all-pass transfer functions.
be a realization of a descriptor system (1) with a regular matrix pencil λE − A, the same number of inputs and outputs, m = p, the system's transfer function G(s) and ς > 0 a real constant. Also, it is assumed that the descriptor system is R-controllable and R-observable. Then G(s) is all-pass, i.e., G(s)G T (−s) = ς 2 I m holds, if and only if the following conditions are satisfied: 1. There are symmetric matrices G pc and G po with 2. The matrices G pc and G po are the solutions of the projected generalized continuous-time Lyapunov equations

The proper Hankel singular values satisfy
4. Let G(s) = G sp (s) + P (s) be decomposed into the strictly proper part G sp and the polynomial 5. Also, the following constraints hold Proof. The proof can be found in the Appendix.
Theorem 1 can be used to derive more general construction formulas for all-pass error systems, see [8] for generalized formulas for an invertible E matrix. Also, it shows that the improper part of the system should not change. For the development of an algorithm, the structure of the reducedorder model (21), obtained from the generalized balanced truncation, can be exploited. So, let the matricesB andĈ be partitioned accordingly to (21) aŝ Using this block partition, the system (Ê,Â,B,Ĉ, D) automatically decouples into its slow subsystemẋ and its fast subsystem First, the fast subsystem (33) is considered. Since the GBT(SR) was used to compute the system (33), there are no zero improper Hankel singular values anymore. As mentioned before and considering Theorem 1, there is no meaningful further reduction concerning the improper Hankel singular values without generating physically meaningless results, so the fast subsystem stays unchanged. Now, let us consider the slow subsystem (32). It is easy to see that (32) is in standard form. Also beneficial properties, resulting from the applied balanced truncation method, still hold for this subsystem which means it is stable and balanced. Let the original system be decomposed as G = G sp + P into its slow subsystem G sp and its fast subsystem P . By the truncation of only zero proper Hankel singular values, the system (32) is a minimal realization of the original slow subsystem G sp . Now, the standard HNA method, mentioned in the previous section, can be applied to (32). As result, an r-th order HNA is computed where E h results from avoiding disadvantageous scaling of the matrices A h and B h . More general transformation formulas for invertible E matrices have been developed in [8]. To get an optimal HNA of the descriptor system (1), now the computed HNA (34) and the reduced-order fast subsystem (33) are coupled In the following theorem, the properties of the resulting GHNA are summarized.
Theorem 2. Let G be a c-stable descriptor system (1) with a regular matrix pencil. The ℓ-th order generalized Hankel-norm approximation (35), with its transfer functionĜ and ℓ = r + ℓ ∞ , has the following properties: 1. The realization ofĜ is conditionally minimal and c-stable.
2. The absolute error in the Hankel semi-norm is given by where ς r+1 (G) is the (r + 1)-st proper Hankel singular value of G.
3. The absolute error in the H ∞ -norm can be bounded by where ς k (G) is the k-th proper Hankel singular values of G.
Proof. Let G = G sp + P be the original system andG = G b + P b the balanced, conditionally minimal realization obtained by the GBT(SR) method. Here G sp , G b denote the slow subsystems and P , P b the fast ones. The GHNA is constructed bŷ where G h is the r-th order HNA (34) of the standard system G b . First, consider part 1. The balanced realizationG is conditionally minimal and c-stable. So by construction (36), both of these properties are transferred to the GHNA. Now consider the error formulas in 2. and 3. Therefor, let E = G −Ĝ be the error system of the GHNA. Then it holds since the balanced realizationG is conditionally minimal and therefor, G b = G sp and P b = P . Using the error bound of the standard method (15) one obtains Using the same approach, the error in the H ∞ -norm is given by if the H ∞ -norm error bound for the standard r-th order HNA from [2] is used.
In Algorithm 1, the complete GHNA method is summarized.

Approximate GHNA
The GHNA method can quickly become numerically unstable. This problem arises from the transformation formulas (12) for the construction of a scaled all-pass error transfer function. It is easy to see that the diagonal matrix Γ =Σ 2 − ς 2 r+1 I nmin−k can lead to large numerical errors for small proper Hankel singular values in further computations. This happens if either the chosen value ς r+1 or the remaining proper Hankel singular values inΣ are very small. One preventive measure was the usage of the descriptor system structure (34) to avoid unnecessary scaling by Γ. In further considerations, only the case of too small remaining Hankel singular values is treated.
Small proper Hankel singular values can arise from numerical errors during the computation of the minimal realization. Therefor, one approach to solve this problem is to compute a smaller balanced truncation of the slow subsystem than the minimal realization such that too small proper Hankel singular values are cut off. In this case, an additional error is made since the balanced realization is only an approximation of the original system. To get a measure for the additional  (6) and (7) for the Cholesky factorizations G pc = R p R T p and G po = L p L T p . 2: Solve the discrete-time Lyapunov equations (8) and (9) for the Cholesky factorizations G ic = R i R T i and G io = L i L T i . 3: Compute the two skinny singular value decompositions 4: Compute the transformation matrices

5:
Compute the minimal balanced realization of the slow subsystem 6: Choose the proper Hankel singular value ς r+1 . 7: Permute and partition the Gramians of the slow subsystem G pc =Ǧ po = diag(Σ, ς r+1 I k ), and the corresponding system matriceš 8: Compute the all-pass transformatioñ where F is anti-stable and G h stable with the realization (E h , A h , B h , C h , D h ). 10: Compute the balanced realization of the fast subsystem 11: Couple the resulting subsystems error, let G b be the computed balanced truncation of order n b of the slow subsystem G sp . Then it has been shown in [11] that in the Hankel semi-norm it holds with n f the order of the slow subsystem G sp . For the overall error, let G = G sp + P be the original descriptor system andG = G b + P b the balanced realization with G b of order n b . The generalized Hankel-norm approximation is denoted byĜ = G h + P b , where the r-th order standard Hankel-norm approximation G h was computed from the balanced realization G b . Using (37) one obtains Concerning the H ∞ -norm, the approach in (38) can be used to get which is the same error bound as for the exact method. This approximate version of the GHNA takes advantage of the use of the GBT(SR) method in form of the adaptive choice of the order n b . It is possible to choose the order n b with respect to the proper Hankel singular value ς r+1 such that In this case, the resulting additional error becomes negligible small concerning the original Hankel semi-norm error. But the corresponding matrix Γ leads to a better conditioned problem. The algorithmic adjustments in the implementation of the GHNA method are small, since only the truncation of non-zero proper Hankel singular values has to be allowed in the generalized balanced truncation method. In this case, the Σ 2 term in (17) with the undesired proper Hankel singular values is not zero and only the matrices U 1 , Σ 1 , and V 1 are used for further computations.
Another advantage of the approximate algorithm can be found in the computation of the balanced truncation. The GBT(SR) method needs to scale the transformation matrices (19) using the inverse remaining Hankel singular values which is more accurate if the small proper Hankel singular values are truncated. Also in the sense of computational costs, this approximate method has advantages. The further steps of the algorithm, i.e., the all-pass transformation and additive decomposition, are extremely costly for large matrices in terms of computational time and memory usage. Therefor, it is advantageous to already have a small balanced realization for the further computations.

Application to Sparse Systems
A frequently appearing case in practice is the model reduction of large-scale sparse descriptor systems. In this case, the system matrices E and A from the descriptor system (1) are in a largescale sparse form, i.e., the dimension n is large, the matrices can be stored using O(n) memory and the matrix-vector multiplication can be computed in O(n) effort. Often such matrices result from the discretization of partial differential equations.
The transformation into a balanced realization does not preserve the sparsity of the system matrices. Therefor, the GHNA method can only be adapted to sparse systems in the first two steps. This concerns the computation of the solutions of the generalized projected Lyapunov equations (6)- (9). It has been observed that the eigenvalues of the symmetric positive semidefinite solutions of Lyapunov equations with low-rank right-hand sides generally decay rapidly. The same result holds for the generalized projected Lyapunov equations [21]. Therefor, the system Gramians can be approximated by low-rank Cholesky factorizations, e.g., G pc ≈ Z pc Z T pc with Z pc ∈ R n×k and k ≪ n.
For the proper system Gramians, the computation is done by adapting existing low-rank methods, e.g., Krylov subspace methods or low-rank ADI methods. In this case, the right-hand side has to be replaced by the projected form from the Lyapunov equations (6), (7). Additionally, it is recommended to project the solution back into the corresponding subspace after some steps of the methods due to a drift-off effect.
In contrast to this, for the improper system Gramians full-rank factorizations can be constructed explicitly such that [21] for more details. Thereby, the size of the full-rank factorizations is bounded by the number of inputs m or outputs p times the system's index ν. This corresponds to the overall bound of the non-zero improper Hankel singular values (16).
Still for using these methods, the spectral projections P ℓ , P r , Q ℓ and Q r have to be computed. But for many problems, these spectral projections can be applied by exploiting the special structure of the problem; see [21] for some examples.

The Projection-Free Approach
In case of unstructured problems, there are no explicit construction formulas for the spectral projectors P ℓ , P r , Q ℓ and Q r , so they have to be explicitly computed for the use in the generalized projected Lyapunov equations (6)- (9). But as for the GBT(SR) method, an alternative approach to the use of spectral projectors can be given; see [19].
As already used in the GHNA algorithm, the GBT method can be interpreted as a decoupling of the original system into the slow and fast subsystems and the individual reduction of both. Therefor, consider the following generalized block triangular form. There are orthogonal matrices U, V ∈ R n×n such that where the matrix pencil λE f − A f contains all the finite eigenvalues of λE − A and the matrix pencil λE ∞ − A ∞ has only infinite eigenvalues. For the computation of a block diagonalization of the system, the coupled Sylvester equations have to be solved for Y and Z; see [5]. Using all of these matrices for the restricted system equivalence transformation where the remaining matrices are constructed as Obviously, the realization in (41) decouples into the fast and slow subsystems of (1). Since the spectral projectors of the subsystems are identity matrices, the corresponding Lyapunov equations (6)- (9) simplify to for the slow subsystem and for the fast subsystem. These Lyapunov equations can be computed without the spectral projections. The matrices X pc and X po correspond to the parts of the proper controllability and observability Gramians, which contain the potentially non-zero proper Hankel singular values. The same holds for X ic , X io and the improper system Gramians. For the rest of the algorithm, only the transformations have to be restricted to the subsystems. The projection-free approach is implemented in the version 3.0 of the MORLAB toolbox; see [7]. In this special implementation, the block diagonalization of the system is done by using a block transformation approach based on the following generalization of Theorem 4.1 from [13].
Theorem 3. Let Γ ⊂ C be a region in the complex plane which contains n 1 eigenvalues of the matrix pencil λE − A. Let Q, Z ∈ R n×n be orthogonal matrices that transform the matrix pencil λE − A into the upper block triangular form: 22 ) = ∅. Similarly, let U, V ∈ R n×n be orthogonal matrices that transform the matrix pencil λE − A into the upper block triangular form: 22 , E Proof. The proof can be found in [22,Section 5.2].
In contrast to the approach above, it is not necessary to compute the solution of the coupled Sylvester equations and, due to the block orthogonal structure of the transformation matrices, the right-hand sides are usually better conditioned than (42). In MORLAB, the right matrix pencil disk function method is used to generate the block transformation matrices, see [22] for more details on the implementation. Additionally, Theorem 3 can be used to compute the additive decomposition in step 9 of Algorithm 1 by separating the eigenvalues with negative and positive real-parts.

Numerical Examples
Two examples have been chosen to demonstrate the introduced GHNA method. All the computations were done on a machine with one Intel(R) Core(TM) i7-6700 CPU processor running at 3.40GHz and equipped with 8 GB total main memory. The computer is running on Ubuntu 16.04.4 LTS and uses MATLAB 9.1.0.441655 (R2016b).

Semi-Discretized Stokes Equation
First, the method is tested on a large-scale sparse example. The Stokes equation describes the flow of fluids at very low velocities without convection and coincides with the linearization of the Navier-Stokes equation around the zero-state. The spatial discretization of the Stokes equation by the finite volume method leads to a descriptor system of the forṁ where v h and p h are the semi-discretized vectors of velocity and pressure, respectively, and the matrices B 1 , B 2 , C 1 , C 2 are all vectors. For matrix pencils like in (43) the spectral projectors P ℓ and P r are given by explicit construction formulas is the orthogonal projector onto the kernel of A T 12 along the image of A 12 ; see [20]. The generation of data is based on the test example 3.3 in [16]. The Stokes equation was discretized on a uniform staggered grid of 80 × 80 points which leads to a descriptor system of the size n = 19, 039, where the matrix pencil λE − A has n f = 6, 241 finite and n ∞ = 12, 798 infinite eigenvalues. The data was generated to get a full-rank A 12 such that the system (43) is of index 2.
For the computation, the implementation of the GHNA method was adjusted to the sparse system case, as described in Section 4.2, and for the solution of the projected continuous-time Lyapunov equations (6) and (7), the solvers from version 1.0.1 of the M-M.E.S.S. toolbox have been used [15]. See the demo file bt mor DAE2.m in [15] for the used parameter settings. With these adjusted solvers, the two iterations for the low-rank factors quickly converged after 31 and 32 iteration steps as shown in Figure 1. An approximation of the non-zero proper Hankel singular values has been computed and plotted in Figure 2 using the low-rank factorizations of the proper system Gramians.
As mentioned before, it is numerically more stable using a balanced truncation of the slow subsystem than the minimal realization. For this reason, a tolerance for the allowed proper Hankel singular values was computed as log(n) · ǫ and multiplied with the largest proper Hankel singular value, with n the order of the system and ǫ the machine epsilon. The resulting bound is also shown in Figure 2 and the computed balanced realization is of order 21.
To compute a fourth order standard Hankel-norm approximation of the slow subsystem, the fifth proper Hankel singular value ς 5 = 1.8370 · 10 −6 was chosen. The additive decomposition of the transformed realization (12) was made by using the ml adtf dss routine from version 3.0 of the MORLAB toolbox [7]. The projected generalized discrete-time Lyapunov equations (8) and (9) were constructed as shown in Section 4.2. In contrast to the continuous-time case, every iteration step was reprojected since the iteration converges after 2 steps at maximum. As result only one non-zero improper Hankel singular value θ 1 = 5.3046 · 10 −18 was computed. This implies that the reduced-order system would be of index 1. In this case, the fast subsystem (33) is equivalent to a feed-through term of the form −C ∞ B ∞ = −1.875 · 10 −17 . Since this value is negligible small compared to the resulting feed-through termD = ς 5 from the GHNA method, the state corresponding to this improper Hankel singular value was truncated, too.
To sum up, the original semi-discretized Stokes equation is approximated by a GHNA of order 4 (r = 4, ℓ ∞ = 0). The error of the original and reduced-order transfer functions in the spectral norm can be seen in Figure 3. Additionally, the corresponding H ∞ error bound as well as a reduced-order model of the same order computed by the GBT(SR) method are plotted to get an impression of the approximation behavior. The shown error behavior of the GHNA is typical. Since the reduced-order model is based on an all-pass error transfer function, the error behavior becomes nearly all-pass if the influence of the anti-stable part is negligible small. Also, the error of the GHNA approaches the chosen proper Hankel singular value ς 5 , which is exactly the error of the approximation in the Hankel semi-norm.
Further examples and tests of the sparse implementation of the GHNA method can be found in [22].

A Damped Mass-Spring System
As a second example, a damped mass-spring system with a holonomic constraint is considered here. The detailed construction of the system can be found in [14]. The vibrations of the resulting system are described by a system of second-order equations where p(t) is the position vector, λ(t) ∈ R is the Lagrange multiplier, K, D ∈ R g×g are the tridiagonal stiffness and damping matrices, M = diag(m 1 , . . . , m g ) is the mass matrix and G = 1, 0, . . . , 0, −1 is the constraint matrix. The input matrix is given by B u = e 1 and three positions of masses are measured by C p = e 1 , e 2 , e g−1 T , where e i is the i-th column of I g .
For the application of the GHNA method, the system (44) has to be rewritten in first-order form. Therefor, the velocity vector v(t) =ṗ(t) is introduced and all states are collected in x(t) = p(t) T , v(t) T , λ(t) T , such that the system (44) can be rewritten in the form This linearization is an index-3 descriptor system. The number of masses was chosen as g = 1500, which leads to n = 3001 states in the linearized system (45). For the computation of the GHNA, the ml hna dss method from version 3.0 of the MORLAB toolbox has been used [7]. In this function, the projection-free approach from Section 4.3 is implemented as mentioned there. For the computation of the additive decompositions, the right matrix pencil disk function is used and the generalized Lyapunov equations are solved via the matrix sign function method; see, for example, [3] and [4]. More details on handling descriptor systems with the MORLAB toolbox can be found in [9]. The computed proper Hankel singular values and the used bound for the minimal realization of the system can be seen in Figure 4.
The computed reduced-order model is of order 6 (r = 6, ℓ ∞ = 0). So also in this case, the reduced-order model does not contain algebraic constraints anymore, which means theÊ matrix is regular. The absolute error of the GHNA is plotted in Figure 5 with the corresponding H ∞ error bound and the error of the GBT(SR) reduced-order model for comparison.

Conclusion
An algebraic characterization of descriptor systems with all-pass transfer function was proven and based on this explanation, an efficient algorithm for the computation of the generalized Hankelnorm approximation was developed by exploiting the generalized balanced truncation square root method. To get a numerically more stable algorithm, an approximate version of the Hankelnorm approximation was introduced. For an efficient practical usage, the introduced method was considered for sparse large-scale systems as well as for unstructured dense systems. The approximation behavior of the method was shown on large-and medium-scale examples.
In contrast to the approach of Cao, Saltik, and Weiland [10], the method, introduced in this paper, has several numerical advantages. It has a more stable and efficient computational behavior, due to the fact that the Weierstrass canonical form does not have to be computed. Also, the introduced method can be applied to more general descriptor systems since C-controllability and C-observability were not assumed.
Both cases are a contradiction to the condition that the coefficients are real and non-zero. Therefor, an all-pass transfer function cannot be improper.
Strictly Proper Case: Now, let's assume that G is a strictly proper transfer function. Using the same argumentation as in the improper case, we get that the product of a strictly proper transfer function with its para-Hermitian is also strictly proper. Now, Theorem 1 can be proven.
Proof. At first, we can assume w.l.o.g. that ς = 1, since the system can be scaled to that case bỹ So the expressions (28) and (29) hold. Since the matrix pencil λE − A is assumed to be regular, there are non-singular matrices Q, Z ∈ R n×n , which transform the matrix pencil into the following block diagonal structure where λE f − A f contains all the finite eigenvalues of λE − A and λE ∞ − A ∞ contains only infinite eigenvalues. These transformation matrices can be used on the complete system as a restricted system equivalence transformation (QEZ, QAZ, QB, CZ, D) This system decouples into its slow and fast subsystem The slow subsystem corresponds to the strictly proper part of the transfer function and the fast subsystem to the polynomial part. Then, the constant part of the transfer function is then given by The equality M −1 0 = M T 0 was already proven above. From the R-controllability and R-observability assumption together with the regularity of E f , it follows that there exist invertible matrices T, W ∈ R n×n which transform the realization of the inverse transfer function into the realization of the para-Hermitian one, with Now, these expressions can be reformulated. From (49) we obtain From (50) we get The equation (48) can be reformulated as And as last one, for (47) it holds Therefor, T and W −T as well as T −1 and W T satisfy the same set of equations, which means that W = T −T . Using this, the expressions (47)-(50) are equivalent to The expressions (51), (52) and (54) give the T as solution of the following system of matrix equations By setting the symmetric matrixG pc = −T E −T f = −E −1 f T T the equation system can be rewritten as Analogously, it follows with the symmetric matrixG po = −T −1 E −1 f = −E −T f T −T . For the matricesG pc andG po , the following matrix product is considered and alsoG As last step in this part, the realization has to be back-transformed into the original one. By multiplying (61) from the left with Q −1 and from the right with Q −T we get