BALANCING BASED MODEL REDUCTION FOR STRUCTURED INDEX-2 UNSTABLE DESCRIPTOR SYSTEMS WITH APPLICATION TO FLOW CONTROL method

. Stabilizing a ﬂow around an unstable equilibrium is a typical problem in ﬂow control. Model-based designed of modern controllers like LQR/LQG or H ∞ compensators is often limited by the large-scale of the discretized ﬂow models. Therefore, model reduction is usually needed before designing such a controller. Here we suggest an approach based on applying balanced truncation for unstable systems to the linearized ﬂow equations usual- ly used for compensator design. For this purpose, we modify the ADI iteration for Lyapunov equations to deal with the index-2 structure of the underlying descriptor system eﬃciently in an implicit way. The resulting algorithm is test-ed for model reduction and control design of a linearized Navier-Stokes system describing von K´arm´an vortex shedding .


1.
Introduction. Stokes, Navier-Stokes or Oseen equations play an important role in describing various problems in fluid dynamics and engineering applications. The spatial discretization of such equations using finite difference or finite element methods produces large-scale structured index-2 descriptor systems [24,5]. In case of controller design, simulation and design optimization, working with these large-scale systems is often problematic due to computational complexity and storage requirements. The idea of model order reduction (MOR) [3,13,28] is to approximate a large scale dynamical system by a substantially lower dimensional system which has nearly the same input-output behavior. This lower dimensional model then, e.g. serves as the basis for feedback control design. The system theoretic method balanced truncation (BT) [3,13] is a particular technique which preserves important properties of the system (such as asymptotic stability) and features an a priori bound on the approximation error.
The theory of BT for large-scale asymptotically stable descriptor systems has been developed first by Stykel in [32,31]. There the author discusses the general framework of the BT method for a descriptor system. In principle, the proposed method is based on splitting the descriptor system into proper and improper subsystems corresponding to the deflating subspaces of the associated matrix pencil with respect to finite and infinite eigenvalues, and on then reducing only the order of the proper subsystem. This approach requires the availability of spectral projectors onto the respective subspaces. Recently, a method for BT of structured large-scale descriptor systems of index 2 has been developed [22] that avoids the computation of spectral projectors. Instead it implicitly performs an index reduction by projection to the inherent or hidden manifold on which the solution evolves.
In this paper, we study the BT based model reduction of an unstable structured index-2 descriptor system. In principle, one can apply the above BT technique to this model by stabilizing the system first using a proper stabilizing feedback matrix (SFM) and then following the approach in [22]. Picking up the ideas in [36] we apply the truncating transformations computed with respect to the stabilized system to the original unstable system. To compute the controllability and observability Gramian factors (which are the main ingredients in determining the balancing transformations), we need to solve two projected algebraic Lyapunov equations of the stabilized system. Again following [36] we employ Bernoulli stabilization to derive the SFM. The main advantage of the Bernoulli stabilization is that it only changes the anti-stable eigenvalues of the system. Thus the required Bernoulli equation can be restricted to these and is of the same dimension as the corresponding eigenspaces.
The main computational bottleneck in BT methods is the efficient solution of large-scale Lyapunov equations, where here, the additional difficulty of the involved (implicit) spectral projection has to be treated. Compared to the presentation in [22], this paper also presents an updated version of the low-rank Cholesky factoralternating direction implicit (LRCF-ADI) algorithm to solve these projected Lyapunov equations. Here, we call this approach the generalized sparse (GS-)LRCF-ADI method. In order to ensure fast convergence of GS-LRCF-ADI, proper selection of the ADI shift parameters is crucial. In Section 6.2 we also discuss and resolve the difficulties in computing shift parameters for the models of flow control considered here.
Stabilization of the non-stationary incompressible Navier-Stokes equations around a stationary solution using a Riccati-based feedback has received considerable attention regarding theory as well as numerical methods during the recent years. In the Riccati-based boundary feedback stabilization technique [6], the most challenging task is to solve the corresponding algebraic Riccati equation (ARE) for the high dimensional model. In this paper we numerically investigate the potential of a feedback control for the original model computed using a reduced-order model. In direct comparison this can be computationally much cheaper if the number of Newton steps for solving the Riccati equation is larger than 2.
The proposed method is applied to data for a spatially FEM semi-discretized linearized Navier-Stokes model. Numerical results are discussed for both the BT model reduction as well as the reduced-order model based stabilization. Figure 1. Initial discretization of the von Kármán vortex street with coordinates, boundary parts and observation points (source [6]).
This article is organized as follows. In Section 2, we describe the model which is considered for our numerical test. The fundamental idea of BT for both stable and unstable systems and the LRCF-ADI methods are repeated briefly in Section 3. Section 4 collocates the ideas for converting the index-2 system into an ODE system applying the projector onto the hidden manifold together with the properties of the projector. In Section 5 we discuss BT based MOR techniques for structured index-2 unstable systems. Section 6 consists of the GS-LRCF-ADI algorithm for the implicit solution of the projected generalized algebraic Lyapunov equations and also the techniques for selecting ADI shift parameters for models of the given structure. In the subsequent section we illustrate numerical results, obtained by applying our approaches to the model introduced in Section 2 with different spatial resolutions and consequently different matrix dimensions.
2. Model example. In his seminal book [30] Sontag, among others, presents the linearization principle. That principle basically states that a general nonlinear model can be stabilized by a linear quadratic regulator (LQR) for a linearization of itself in the vicinity of the linearization point. The basic idea is that if the regulator is working properly, the vicinity where the linearization is a proper approximation of the nonlinear system is never left. This principle has been exploited by the authors in [6] for a Navier-Stokes model for the von Kármán vortex street. The linearized Navier-Stokes equations that arose there and that we consider in this paper are where v, w denote velocity vectors, p the pressure and Re is the Reynolds number. The vector w represents the stationary solution of the incompressible nonlinear Navier-Stokes equations and v is the deviation of the original state from the stationary solution. The boundary and initial conditions, as well as the derivation of this model, are given in [6]. There the authors apply a mixed finite element method based on the well known Taylor-Hood finite elements [23] to discretize equation (1). The coarsest discretization of the domain for the von Kármán vortex street example from [6] is shown in Figure 1. This yields the differential-algebraic equations where v(t) ∈ R n1 denotes the nodal vector of discretized velocity deviations, p(t) ∈ R n2 the discretized pressure, u(t) ∈ R m are the inputs, and E 1 , A 1 ∈ R n1×n1 , A 2 ∈ R n1×n2 , B 1 ∈ R n1×m are all sparse matrices. Additionally, the vertical velocities in the observation nodes depicted in Figure 1 in the domain are modeled by the output equation with the output y(t) ∈ R p and the output matrix C 1 ∈ R p×n1 . The above system remains stable, i.e., the finite spectrum of the matrix pencil is located in the negative half plane C − , as long as the Reynolds number Re is small. However, already for moderate Reynolds numbers (e.g., in the configuration of Figure 1 Re ≥ 300) a few finite eigenvalues move to the positive half plane, C + [2]. The key problem in the LQR approach for the model under investigation is to compute the boundary feedback stabilization matrix K 1 (see e.g., [6]), such that the stabilized system has the following form: The authors in [6] presented an algorithm (see [6,Algorithm 2]) to compute K 1 which is based on the standard linear-quadratic regulator approach [30,16] for a projected semidiscretized state-space system. The most challenging part in this algorithm is to solve the usually very large, generalized, projected algebraic Riccati equation (GARE) based on the full order semidiscretized model. In this paper our aim is to replace this full order system (2) together with its associated output equation (3) by a substantially lower dimensional model where the responses (both in frequency and time domain) of the reduced-order model are good approximations to those of the full model. Note that for balanced truncation, we can always guarantee thatÊ = I, which we thus assume in the following. We employ the reduced-order model to compute an approximation to the optimal LQR feedback matrix of the full system. The main advantage of this approach is that we only need to solve two projected algebraic Lyapunov equations in order to derive the reduced-order model instead of one Lyapunov equation per Newton step in the solver for the GARE, which are usually many more [14].
Based on the reduced model (6) the GARÊ is now much smaller in dimension. It can thus easily be solved forX using classical solvers as, e.g., the MATLAB care command. The stabilizing feedback matrix for the reduced model (6) then isK The ROM-based approximation to the stabilizing feedback matrix for the full order model can now be retrieved as where T L is the left balancing transformation (see Section 5) used to compute the reduced-order model.

3.
Background. In this section we briefly review the basic idea of standard balanced truncation for both stable and unstable systems. We also repeat the basics for the LRCF-ADI method for computing the required Gramian factors in large scale settings, to put our new contributions for the index-2 case in Algorithm 2 (see Section 6) into perspective.
3.1. Standard BT and Low-Rank-ADI. Consider a generalized linear timeinvariant (LTI) continuous time system in which E ∈ R n×n is non-singular, and A ∈ R n×n , B ∈ R n×p , C ∈ R m×n and D ∈ R m×p . Assume that the original system is asymptotically stable, i.e., all eigenvalues of the pair (A, E) lie in C − . Then the controllability Gramian P ∈ R n×n and the observability Gramian Q ∈ R n×n , which are the solutions of the two Lyapunov equations (see, e.g., [3]) respectively, are symmetric and positive semi-definite, and therefore have symmetric decompositions P = RR T and Q = LL T . (10) The balancing transformation can be formed using the SVD and defining Here, U 1 and V 1 are composed of the leading k columns of U and V, respectively, S 1 is the first k × k block of the matrix S = diag (σ 1 , σ 2 , . . . , σ k , . . . σ n ), where σ i , i = 1, 2, . . . , n are the Hankel singular values of the system (9) [3,18]. Now the reduced-order model (ROM) is derived aŝ The ROM obtained in this way satisfies the global error bound [3,18] where G(s) := C(sE − A) −1 B andĜ(s) :=Ĉ(sÊ −Â) −1B are called the transfer functions [3] of the full and reduced model, respectively, and . H ∞ denotes the H ∞ -norm. The relation (14) is indeed an a priori error bound. Hence, one can easily determine the required dimension of the ROM for a given error tolerance.

11:
Update low-rank solution factor In [33] this procedure for the BT method is summarized as the square-root (SR) algorithm.
If the numbers of inputs and outputs are very small compared to the number of degrees of freedom (DoFs), then usually the Gramians P and Q can be expressed in low-rank factored form [26,4,20], i.e., R and L are thin rectangular matrices. They can for example, be computed by the LRCF-ADI iteration [26,25,12]. Recently, the formulation of the LRCF-ADI method has been updated in [9], such that real Gramian factors can be computed by clever handling of complex shift parameters. Moreover the same authors [10] introduced an efficient technique to compute a lowrank factor of the ADI residual as an integral part of each iteration step and thus cheaply evaluate stopping criteria for the algorithm.
Combining the ideas from [9] and [10] the LRCF-ADI is summarized in Algorithm 1. In this Algorithm either (Ě, Note that if i max > J, then the shift parameters are repeated cyclically.

BT for unstable systems.
This section is concerned with BT for unstable systems via Bernoulli stabilization. Consider the case where the system in (9) is unstable, i.e., some eigenvalues of the system lie in C + . The helpful feature of our investigated example is that still the number of anti-stable eigenvalues is very small. This is exactly the property we exploit for fast computation of the ROMs and ROM based feedback matrices. In the previous subsection we have recalled classical (Lyapunov based) balancing for stable systems. The main ingredients there are the two Gramians (e.g., [3]) which obviously do not exist if the system's unstable poles are controllable. This, however, is exactly the desired case in our motivating example (see Section 2). In [36], the authors use the frequency domain representations of these integrals to extend the definition of the Gramians to systems with no poles on the imaginary axis.
Following the theory in [36] the generalized controllability and observability Gramians P s , Q s for such systems can be computed by solving the algebraic Lyapunov equations where K c = B T X c E and K o = EX o C T are called Bernoulli stabilizing feedback matrices, due to the fact that the matrices X c and X o are the respective stabilizing solutions of the generalized algebraic Bernoulli equations Now, since the Bernoulli stabilization only mirrors the anti-stable eigenvalues across the imaginary axis [7], it is sufficient to solve these Bernoulli equations only on the invariant subspaces corresponding to those eigenvalues. That is, for orthogonal matrices V and W spanning the left and right eigenspaces corresponding to the anti-stable eigenvalues, respectively, we define the Petrov-Galerkin projected system (Ě,Ǎ,B,Č) byĚ Then we solveĚ The projected Bernoulli equations in (17) are solved by the Matrix Sign Function method presented in [7].
The low-rank factors of P s and Q s can also be computed by solving (15) using Algorithm 1, but avoiding to form the closed loop matrices stays crucial. We discuss this issue in more detail in Section 6. Then using the projection matrices as in (11) and applying them to the original unstable system we can compute an unstable ROM that satisfies the error bound as in (14), but with the H ∞ -norm replaced by the L ∞ -norm. Therefore, this error bound can not be translated to a global time domain error bound, as in the classic setting, due to the lack of a Parseval-identitylike result. In fact in the numerical experiments we observed that the error may very well grow over time. 4. Reformulation of the structured dynamical systems. The purpose of this section is to convert the index-2 system (2) into an equivalent ODE system in order to make it fit into the framework for BT based model order reduction as discussed in the previous section. The idea is exactly the same as in [22,Section 3]. For convenience we briefly review it here again. Let us consider a projector of the form Clearly one then has These properties imply i.e., the image of Π T is exactly the subspace where the algebraic condition (2b) of the DAE investigated here is satisfied. Now applying the projector in (2) and the output equation (3), and thus reformulating the system on this exact subspace, we obtain the following projected system The system dynamics of (21) are projected onto the n m := n 1 − n 2 dimensional subspace Range (Π T ) [22]. This subspace is, however, still represented in the coordinates of the n 1 dimensional space. The n m dimensional representation can be made explicit by employing the thin singular value decomposition (SVD) where Θ l = U 1 and Θ r = V 1 and in which U 1 , V 1 ∈ R n1×nm consist of the corresponding leading n m columns of U , V ∈ R n1×n1 . Moreover, Θ l , Θ r satisfy This representation is always possible since Π is a projector and therefore, Σ 1 = I nm (n m dimensional identity matrix). Inserting the decomposition of Π from (22) into (21) and consideringṽ = Θ T l v, we get System (24) practically is system (9) with the redundant equations removed by the Θ r projection. We observe that the dynamical systems (2), (21) and (24) are equivalent in the sense that their finite spectra are the same [17, Theorem 2.7.3] and the input-output relations are the same, i.e., they realize the same transfer function.

5.
BT for index-2 unstable descriptor systems. In this section we show how to avoid forming (24) explicitly. Suppose that we want to apply balanced truncation to the system (24). To this end, we need to solve the Lyapunov equations where are the corresponding projected controllability and observability Gramians. Again, K 1 c and K 1 o are the Bernoulli stabilizing feedback matrices and can be computed as described in Section 3.2. The solutionsP ,Q of (25) are unique since we assured that the respective dynamical system is asymptotically stable and symmetric positive (semi-)definite since the right hand side is semi-definite. Now multiplying (25) by Θ l from the left and Θ T l from the right and exploiting that Θ r = Π T Θ r (e.g., due to (22), (23)) we obtain where P = Θ rP Θ T r and Q = Θ rQ Θ T r . These are the respective controllability and observability Lyapunov equations for the realization (21) and the solutions P, Q ∈ R n1×n1 are the corresponding controllability and observability Gramians. The system (26) is singular due to the fact that Π is a projection. Uniqueness of solutions is reestablished by the condition that the solutions satisfy P = Π T P Π and Q = Π T QΠ.
It is also an easy exercise to go back to (25) from (26). Let us consider P ≈ RR T , Q ≈ LL T andP ≈RR T ,Q ≈LL T . Then R, L,R andL are called approximate low-rank Cholesky factors. They fulfill the relation For large-scale problems, however, computing Θ r is usually impossible due to memory limitations. Therefore, R and L are computed by solving (26). The balancing transformations for (24) arẽ where U k , V k , ∈ R nm×k consist of the corresponding leading k columns of U, V ∈ R nm×nm , and S k ∈ R k×k is the upper left k × k block of S in the SVD Observing further that R T ΠE 1 Π T L =R T Θ T r E 1 Θ rL = U SV T , the projection matrices for the system (21) can be formed as As in [22] we find that Now we apply the transformations T R and T L in (21) to find the reduced-order model as in (6) wherê (29) we can avoid the explicit usage of Π and find Eventually, we see that the reduced-order model (6) is obtained without forming the projected system (21). In the next section we will show how to compute R and L using a tailored version of the LRCF-ADI iteration without using Π explicitly.
6. Solution of the projected Lyapunov equations. In order to apply the aforementioned balancing based MOR we need to solve the projected Lyapunov equations (26). We have seen above that the Π T invariant solution factors enable us to compute the corresponding truncating transformations. The approach here is different from the spectral projection based work by Stykel in that we are applying the E 1 -orthogonal projection to the hidden manifold, where Stykel uses the orthogonal projection (in the Euclidean sense) to the eigenspaces corresponding to the finite poles of the system. In fact both methods project to the same subspace only considering orthogonality in different inner products. Here we are concerned with two main issues. First we discuss the reformulation of the basic low-rank ADI Algorithm for the projected Lyapunov equation that ensures the invariance of the solution factor and the computation of the correct corresponding residual factors. We are lifting the ideas of [22] to the reformulation of the LR-ADI iteration in Algorithm 1. For the spectral projection methods the analogue procedure has been discussed in [14]. In the second part we address the important issue of ADI shift parameter computation. There the main issue in the DAE setting is to avoid the subspaces corresponding to infinite eigenvalues in order to correctly compute the large magnitude Ritz values involved in many parameter choices. The crucial point in both parts is to provide methods that use the projection Π T at most implicitly and never form the projected system (21).
6.1. GS-LRCF-ADI for index-2 unstable systems. This section is concerned with the efficient solution of the Lyapunov equations in (26) to compute the lowrank Gramian factors using the LRCF-ADI method as discussed in Section 3. Here, we consider the projected controllability equation elaborately. The observability equation can be handled analogously. For convenience we rewrite the Lyapunov equations (26) asÃPẼ T +ẼPÃ T = −BB T , In the i-th iteration step of the ADI method, the residual of the controllability Lyapunov equation (32) can be written as

11:
Update low-rank solution factor To compute the low-rank controllability Gramian factorR we follow Algorithm 1. In the i-th iteration step, V i is computed from which enables us to update the residual factor according tõ In complete analogy to [22,Lemma 5.2], we observe that instead of solving (33), one can compute V i by solving where the special case i = 1, here especially the computation of the initial residual factorW 0 , is discussed in detail below.
In this equation the inversion of (A − B K) should in practice be computed using the Sherman-Morrison-Woodbury formula (see, e.g. [19]): to avoid explicit forming of the (usually dense) matrix A − B K. Note that also A −1 should never be formed, but implemented as a sequence of linear system solves properly exploiting the sparsity of A. That means, e.g., the term A −1 B is to be understood as the solution of the linear system of equations AZ = B, possibly with several right-hand sides depending on the number of inputs, i.e., columns of B. This can be achieved by any suitable sparse direct or iterative solver. In accordance with [22,Lemma 5.2], again the computed V i in (35) satisfies V i = Π T V i . Therefore, the correct projected residual factor in (34) can be obtained bỹ In order to really compute the correct residual, the initial residual must be computed asW 0 = ΠB 1 to ensureW 0 = ΠW 0 . This can be performed cheaply using the following observation.
Again applying the properties in (20), we have A T 2 Ξ = 0. These two relations give (39). Conversely, we assume (39) holds. From the first block row of (39) we get , and thus from the second row it follows Inserting this in the first block row we get as desired This especially ensures Ξ = Π T Ξ, since and thus usingW 0 = E 1 Ξ, we get the desired invarianceW 0 = ΠW 0 . The above findings on the residual factor are summarized in the following lemma.
Lemma 6.2. The residual factor in every step of Algorithm 2 fulfills the relatioñ The whole procedure to compute the low-rank factor of the controllability Grami-anR is summarized in Algorithm 2. Analogous to the derivation in [22], our algorithm computes the correct solution factor. In contrast to the version presented there, we guarantee to compute a real solution factor even if the shifts occur in complex conjugate pairs and we have the low rank residual factors in hand to evaluate stopping criteria cheaply. Still one issue remains open that has not been tackled in the original paper [22]. The shifts that guarantee fast convergence of the algorithm are closely related to the spectrum of the original pencil. The question how these can be computed is answered in the next section. 6.2. ADI-parameter selection. The appropriate shift parameter selection is one of the crucial tasks for fast convergence of the GS-LRCF-ADI iteration. Recently, most of the papers followed the heuristic procedure introduced by Penzl [26] to compute sub-optimal ADI shift parameters µ i , i = 1, 2, . . . , J, for a large-scale problem. Very recently, new shift computation ideas considering adaptive and automatic computation of shifts during the iteration [11,35] have come up. While Penzl's method aims at minimizing the spectral radius of the iteration matrix, restricting the optimization problem to a discrete subset of the left complex half plane, rather than the entire half plane, the new method tries to approximate the optimality of certain projected eigenvalues [8] in the Riemannian manifold optimization framework [34] using a sub-optimal projection basis. The main advantage of the second method is that it is entirely automatic with guaranteed convergence for a large class of systems, while the Penzl method needs some experience in the choice of the discrete set to achieve good results. We present the basic ideas to adapt both the classic and the new methods to our framework in the following two paragraphs.
6.2.1. Heuristic shift selection. The main ingredient of the heuristic method is the computation of the discrete subset of the left complex plane by a number of large and small magnitude Ritz values approximating the system poles. In the case of DAE systems, the computation of Ritz values of large magnitude is causing difficulties due to the existence of infinite eigenvalues. We need to make sure that the infinite eigenvalues are avoided. This can be achieved by the following fact that is a direct consequence of [15, Theorem 3.1]. Corollary 6.3. The matrix pencil transforms all infinite eigenvalues of the pencil P(λ) (see (4)) to 1 δ while at the same time preserving its finite eigenvalues.
Thus from the pencil P δ we can compute the desired Ritz values of large magnitude via an Arnoldi iteration [27]. The parameter δ can easily be determined after the small Ritz values β i have been computed with respect to the original pencil. In order to ensure that 1 δ is avoided by the Arnoldi process for the large magnitude Ritz values, and thus only finite eigenvalues of the original pencil are approximated, one could, e.g., set δ = 1 min i Re (βi) . For the unstable case the corollary obviously has to be applied with A 1 replaced by A c .

Adaptive shift selection.
A second shift computation strategy we use in the numerical experiments follows the lines of the adaptive shift strategy proposed in [11]. There, the shifts are initialized by the eigenvalues of the pencil projected to the span of W 0 . Then, whenever all shifts in the current set have been used, the pencil is projected, e.g., to the span of the current V i and the eigenvalues are used as the next set of shifts. Here, we use the same initialization. For the update step, however, we extend the subspace to all the V i generated with the current set of shifts and then choose the next shifts following Penzl's heuristic with the Ritz values replaced by the projected eigenvalues computed with respect to this larger subspace. Note that in lack of the conditions for Bendixon's theorem, we cannot guarantee that the projected eigenvalues will be contained in C − . Should any of them end up in the wrong half-plane, we neglect them. Should the resulting set  Table 2. The performances of the heuristic and adaptive shifts in the GS-LRCF-ADI iteration.
of shifts become empty due to this, we reuse the previous set as in [11]. Note further that the problem with the infinite eigenvalues does not exist in this case.
Since we have Π T Z = Z, for any orthogonal bases U of a set of columns of Z, we also have Π T U = U . As an immediate result of this fact, the projected pencil (U T A 1 U − λU T E 1 U ) automatically resides on the hidden manifold and can thus only have finite eigenvalues.
7. Numerical results. To assess the performance of the techniques, this section discusses some numerical tests. The method is applied to a set of linearized semidiscretized Navier-Stokes equations as described in Section 2. All the computations were carried out using MATLAB R 7.11.0 (R2010b) on a board with 2 Intel R Xeon R X5650 CPUs with a 2.67-GHz clock speed. The authors of [6] generate different sized models using Reynolds number Re = 500. Table 1 shows the different sizes of the model and distinguishes the dimensions n 1 of the velocity vector (differential variable) and n 2 of the pressure vector (algebraic variable). In all the sets B 1 ∈ R n1×2 and C 1 ∈ R 7×n1 . For Reynolds numbers of 400 and more the described linearized model is unstable. Thus, especially the Reynolds number 500 case discussed here is unstable. The Bernoulli stabilizing feedback matrices K 1 c and K 1 o for all models are computed applying the procedure from [2] and [6,Section 2]. It uses 2 calls of the MATLAB eigs implementation of the Arnoldi method (for P δ and P T δ )to compute the rightmost eigenvalues together with their left and right eigenvectors.
We apply the GS-LRCF-ADI iteration (Algorithm 2) to all aforementioned models to compute the low-rank factorsR andL considering the tolerance 10 −6 . We investigate the performances of both the heuristic and adaptive shifts to implement  Table 4. Balanced truncation tolerances and dimensions of reduced models.
this algorithm. The results are shown in Table 2. For all models we chose 30 optimal heuristic shifts out of 10 large and 80 small magnitude Ritz-values. In the case of the adaptive shifts, in each cycle, 10 proper shift parameters are selected following the procedure discussed above. For computing the initial shifts, first we project the pencil (A 1 − λE 1 ), onto the column space of an n 1 × 100 random matrix. For all the models, looking at the iteration numbers and the execution times, the performance of the adaptive shifts is much better than that of the heuristic shifts. The performances of the heuristic and adaptive shifts are also depicted in Figure 2 for the largest dimensional system Mod-5. This figure shows the convergence rate of the low-rank controllability and observability Gramian factors with respect to iterations (Figures 2a, 2b) and time (Figures 2c, 2d). In both cases the convergence rate for the adaptive shifts is much faster than the heuristic shifts. Note that we use the Frobenius norm to compute the residual norm. We compute the reduced-order models for all the data sets. The dimensions of the original and reduced models are shown in Table 3. If nothing else is stated, the truncation tolerance is set to 10 −5 . The dimension of the ROM can, however, be decreased by increasing the tolerance if desired or required. This is shown in Table 4. Since the numerical results are all comparable, we exemplary present only selected plots (e.g., for the largest model Mod-5).
7.1. Model reduction of the unstable system. This section reviews the numerical experiments for the unstable case. Here we present both frequency and time domain error analyses. The frequency domain error analysis is shown in Figure 3.
In Figure 3a we see the frequency responses of the full and 184 dimensional reducedorder models for Mod-5 with a nice match in the eyeball norm. The absolute and relative deviations between full and reduced-order models are shown, respectively, in Figures 3b and 3c. Here, we can see that the absolute error is bounded by the prescribed truncation tolerance of 10 −5 . For higher frequencies the relative error is  slightly increasing since the frequency response is decreasing more rapidly than the absolute error can. Figure 4 depicts time domain simulation of full and reduced-order models for Mod-5. This figure shows the step responses from Input 1 to Output 1 together with their absolute deviations. To compute the step response we use an implicit Euler method with fixed time step size 10 −2 . Initially, the imposed control is kept inactive, therefore the responses for both (full and reduced) models are constant within the range 0 to 15s. Switching the control to constant unit actuation on Input 1, the responses are oscillating with increasing amplitude in the higher time domain caused by the instability of the model. Here we also see the issue with the balanced truncation error bound for unstable systems since the absolute error is increasing gradually with increasing time.

7.2.
Numerical Experiments for the stabilized system. In Section 2 we have mentioned that the stabilizing feedback matrix for the full model can be computed from the reduced-order model. To this end, we solve the corresponding algebraic Riccati equation (7) arising from the linear-quadratic regulator approach using the MATLAB care command and compute the optimal stabilizing feedback matrix K (see, e.g., [16,Chap. 10]) for the reduced-order model (6). The ROM based approximation to the stabilizing feedback matrix for the full order model (5) is then given by K 1 = K r T T L E 1 . Figure 5 shows the step response (from 1st input to    Step responses of 1st input to 1st output of full and reduced-order models and respective absolute deviations. 1st output) of closed loop full and reduced-order models and their absolute error. For the generation of the step response the same procedure has been followed as for the unstable case above. Note that for a stabilizing feedback the step response  Step responses of 1st input to 1st output of stabilized full and reduced-order models and respective absolute deviations.
system has to be viewed as that of an asymptotically stable system with a constant source term. Thus the output settles at constant non zero values.

Conclusions and Outlook.
We have discussed the model reduction problem for unstable descriptor systems arising from flow control problems. In particular, we have considered linearized Navier-Stokes equations. Their spatial discretization by finite elements leads to index-2 DAEs. This causes some technical difficulties for the application of model reduction based on balanced truncation. Following [22], we have treated the problem by using projectors onto deflating subspaces corresponding to finite eigenvalues of the associated matrix pencils. The explicit application of the projectors is avoided by a re-formulation of the ADI method used to compute approximate low-rank factors of the system Gramians which is the bottleneck in the numerical implementation of the balanced truncation method. Here we show how this can be done for unstable systems, and we also incorporate a number of recent improvements into the resulting ADI method. As an illustrative example, we apply our algorithms to the linearization of the von Kármán vortex shedding at a moderate Reynolds number. It is demonstrated that the resulting reduced-order model can be used to accurately simulate the unstable linearized model and to design a stabilizing controller. Future work will include the realization of the resulting control law for the full nonlinear model. Note that although we restricted our presentation to ADI based balanced truncation, the ideas presented here can be extended to rational Krylov subspace projection based solvers for the Lyapunov equations, in a straight forward way following the lines of [1,21]. Although both treat the iterative rational Krylov Algorithm for model reduction, the rational Krylov spaces can be formed accordingly. Solvers like the extended Krylov subspace method [29] for the Lyapunov equation may require the incorporation of the techniques we presented in the discussion for the Ritz value computation, together with the transformation of the eigenvectors originally presented in [15]. The fair comparison of appropriately sophisticated implementations of all these methods is a pending task, as well.