This paper introduces and analyzes an efficient algorithm for the solution of the generalized continuous algebraic Riccati equation (GCARE) associated with the solution of linear quadratic regulator (LQR) problems governed by Hessenberg index-2 Differential Algebraic Equations (DAEs). This problem arises, e.g., in the computation of feedback laws that stabilize Navier-Stokes flows. The numerical solution of the Riccati equation is challenging because the underlying systems are large-scale and because of the presence of algebraic constraints in the DAE system. To overcome these challenges we extend our inexact low-rank Newton-ADI method developed in [10] for problems governed by ordinary differential equations (ODEs) to this DAE case. The main idea is to use the structure of the Hessenberg index-2 DAE and apply the discrete version of the Leray projector (see Heinkenschloss et al. [16] and Bänsch, et al. [4]) to transform the LQR problem governed by the DAE into a classical LQR problem governed by an ODE. In principle, the standard LQR and Riccati theory as well as the inexact low-rank Newton-ADI method developed in our previous paper [10] can be applied to this ODE problem. This, however leads to a solution approach that is not practical because the projected systems are large-scale and, because of the projection, dense. To arrive at an efficient algorithm, the computations must be presented in terms of the original large-scale sparse system and the structure of the governing DAE system must be exploited. This is done in this paper. In addition, numerical experience with our new algorithm indicates that our control of inexactness and the line-search leads to a start-up phase that reaches the quadratic convergence region of the Newton iteration faster and tends to avoid marginally stable subproblems during intermediate iterations.
The LQR problem and associated Riccati equation considered in this paper have also been solved by Bänsch et al. [4]. However, the focus of [4] was the computation of feedback laws for Navier-Stokes flows, and a basic version of an inexact low-rank Newton-ADI method was applied. Our paper focusses on the solution of the Riccati equation and incorporates many recent improvements. As a result, the algorithm in this paper delivers an approximately 90-times speed-up over the algorithm used in [4]. Benner and Stykel [6] study the solution of projected Riccati equations, which are associated with DAEs. They use so-called spectral projectors, which project onto the right and left deflating subspaces. While these projectors can be applied to general DAEs defined by a regular pencil, in the general case "the projectors [. . . ] are required in explicit form [and the] computation of these projectors is, in general, very expensive" [6, p. 590]. The projector used in our paper is specially designed for the index-2 DAE system arising for fluid flow problems and our Kleinman-Newton-ADI method contains many improvements not yet available in [6]. In principle it is possible to use rational Krylov subspace projection methods (see Simoncini et al. [25,26]) to solve the Riccati equations, but extensions of this approach to the DAE case and numerical comparisons of the latest versions of both approaches are not yet available.
As pointed out above, a main ingredient for the efficiency of our approach is the exploitation of the special structure of the Hessenberg index-2 DAE, in what is called implicit index-reduction. Specifically, we can use structured projectors, rather than generic and expensive spectral projectors. Implicit index-reduction can also be applied to other structured DAE systems, see e.g. [11,14,15,24]. We demonstrate our approach on a large-scale Riccati equation associated with the stabilization of Navier-Stokes flow, but the extension of the techniques described in this paper to other saddle point structured DAEs is straightforward.
This paper is organized as follows. The next section, Section 2, introduces the LQR problem, uses projection onto the constraint manifold to derive a projected Riccati equation, and reviews existence results for both the projected Riccati equation and the LQR problem. Section 3 reviews the main components of our algorithm in [10] applied to the projected GCARE and Section 4 carefully exploits the special structure of the projected GCARE for an efficient numerical realization of the inexact low-rank Newton-ADI method. Finally, Section 5 illustrates the performance of our algorithm on a large-scale Riccati equation associated with the stabilization of Navier-Stokes flow around a cylinder -a problem also solved by Bänsch et al. [4]. As mentioned earlier, the algorithmic improvements in this paper lead to approximately 90-times speed-up over the algorithm used in [4].
Notation. Throughout the paper we consider the Hilbert space of matrices in R n×n endowed with the inner product M, N = tr M T N = n i, j=1 M ij N ij and the corresponding (Frobenius) norm
1/2 . Furthermore, given real symmetric matrices M, N, we write M N if and only if M -N is positive semi-definite, and M N if and only if M -N is positive definite. The spectrum of a symmetric matrix M is denoted by σ (M).
In this section we present the mathematical statement of the LQR problem and the governing Hessenberg index-2 DAE, and we show how it can be transformed into a 'standard' LQR problem governed by an ODE using a projection onto the constraint manifold of the original DAE. Then we apply classical LQR theory to this transformed problem to compute, under standard conditions on the system, the solution of the LQR problem via the GCARE. As mentioned before, the problem transformation is performed to derive the solution, but the computations are done using the original DAE framework. The projection used to convert the DAE into an ODE was first used in a different context by Heinkenschloss et al. [16]. For DAEs derived from a finite element discretization of the Stokes or linearized Navier-Stokes system, Bänsch et al. [4] show that this projection is a discrete version of the Leray projector. Projections have also been used by Benner,Stykel [6] to formulate and solve GCAREs associated with index-2 DAEs, although, as already noted in the introduction, the projection there is different. Except for some extensions in problem statement and notation the material in this section is mostly known from [4,16], but is needed to provide the necessary background that allows us to switch between expressions using the original DAE system and the corresponding expressions using the transformed ODE system. Compared to [4], this section also provides a more detailed link between the representations of the optimal control of the LQR problem derived using the original DAE and transformed ODE system.
where for given u ∈ L 2 (0, ∞), the function y ∈ L 2 (0, ∞) is obtained as the output of the Hessenberg index-2 Differential
Algebratic Equation system
)
(2.2c)
To ensure well-posedness of the LQR, we will make additional assumptions on the system (2.2) in Section 2.2. In the cost functional, we may replace the Euclidian norms by any weighted norm induced by positive definite matrices Q y and Q u . Here, we set both weighting matrices to the appropriate identity for ease of notation. It is straightforward to include non-identify weighting matrices into the problem description and the computational framework.
The LQR problem (2.1), (2.2) arises, e.g., in feedback stabilization of the Navier-Stokes equations, see Bänsch et al. [4] or Raymond [23]. In this context, (2.2a), (2.2b) correspond to the linearized discretized Navier-Stokes equations, and v, p correspond to velocity and pressure, respectively. The problem also arises in feedback stabilization of multi-field flow problems, see Bänsch et al. [3]. In this case, (2.2a) includes additional equations such as linearized reaction equations, and v corresponds to velocities and the other fields, such as concentrations.
If we define
(2.4b)
The structure of (2.2) can be used to convert the LQR problem (2.1), (2.2) into a classical one governed by an ODE. We proceed as in [4,16]. The constraint (2.2b) and the variable p can be eliminated from (2.2a), (2.2b) via the projection
(2.5)
The matrix Π obeys Π 2 = Π and Π M = MΠ T , i.e., it is in fact an M-orthogonal projection. Furthermore, null(Π T ) = range(M -1 G) and range(Π T ) = null(G T ), (2.6) which means that 0 = G T v(t) if and only if v(t) = Π T v(t).
We use the latter property to enforce (2.2b) and multiply (2.2a) by Π to arrive at
(2.7b) If needed, the function p can be computed from v, u using
(2.8) Equation (2.8) is obtained multiplying (2.2a) by G T M -1 and using (2.2b).
Since Π MΠ T ∈ R n v ×n v has an n p -dimensional null-space and cannot be inverted, (2.7) is still not an ODE. However, Π T v(t) ∈ R n v is contained in the n vn p dimensional subspace range(Π T ) and we can explicitly express Π T v(t) as an element of this subspace. This is done using the decomposition
with l , r ∈ R n v ×(n v -n p ) . In particular range( r ) = range(Π T ).
(2.10)
(2.11)
Using the decomposition (2.9), we define
and write the descriptor system (2.7) as
(2.13a)
(2.13b)
The DAE system (2.2) is equivalent to system (2.13), which is an ODE system since with M being symmetric and positive definite so is M, by x T Mx = ( r x) T M( r x) > 0. Furthermore, the LQR problem (2.1), (2.2) is equivalent to the classical LQR problem (2.1), (2.13). We summarize this result in the following proposition.
where
is the unique stabilizing solution of the GCARE
(2.15)
See, e.g., Lancaster, Rodman [19].
The unique stabilizing solution of the GCARE is obtained by applying Newton's method to find a root of the quadratic operator
(2.16)
Given an approximate root X (k) , the new approximation is computed as the solution of
(2.17)
This method is known as the Kleinman-Newton method. See the original paper by Kleinman [17] or the book by Lancaster, Rodman [19].
The system (2.17) is a Lyapunov equation and for large-scale problems the exact Kleinman-Newton method which is defined by (2.17) is impractical. This is particularly true for the Riccati equation (2.15) which is obtained from a large-scale DAE by projection. The projected matrices in (2.12) are not only large-scale, but because of the projections they are also dense. To overcome these difficulties, we need to 'undo' the projections in the numerical computations. We will discuss the details of our solution approach in the next section. In the remainder of this section we provide basic relationships between quantities for the projected problem and quantities for the original problem.
The Kleinman-Newton method applied to the projected Riccati equation (2.15) generates iterates 0
and corresponding feedback matrices
(2.18)
We want to write the corresponding feedback law -(
and
The convergence of the (exact) Kleinman-Newton method can now be expressed in the unprojected variables and in the context of the (2.2). First we show that the stability (detectability) of the system (2.13) is equivalent to the stability (detectability) of the system (2.2).
1. A matrix pencil (A, M) is called stable if it is regular and all the finite eigenvalues of (A, M) lie in the open left half-plane.
2. Let A, B, M be given by (2.12). The triple (A, B; M) is stabilizable if there exists a matrix K ∈ R n v ×n u such that all finite eigenvalues of the matrix pencil
With these preparations, the following result is an immediate consequence of the classical Kleinman-Newton convergence result [17], [19]. See [27,Thm. 4.5] for a detailed proof.
and there is a constant κ such that X (k+1) -X ( * )
r is the solution of the Riccati equation specified in Theorem 4 and
In principle, the large-scale projected GCARE (2.15) can be solved using the Kleinman-Newton method [17]. However, the size and special structure of (2.15) require the inexact solution of the Newton equation, a Lyapunov equation, in each step of the Kleinman-Newton method. Moreover, the explicit use of the large, dense projected matrices M, A, B, C (2.12) must be avoided in computations and the final algorithm must operate with the sparse matrices M, A, B, C (2.3) instead. To adopt our approach from [10] to efficiently solve the large-scale projected GCARE (2.15), we first need to review the main components of our approach in [10].
Our approach in [10] is based on an inexact Kleinman-Newton method with line search. Although the exact and, under additional conditions, inexact Kleinman-Newton method converges with step size fixed to one (see, e.g., Kleinman [17] or Feitzinger et al. [13]), variable step sizes can hugely improve the performance (Benner, Byers [5], Benner et al. [10]). We will also observe this in our numerical tests, see Fig. 2 in Section 5. The line search method and analysis in [5] are based on exact Lyapunov equation solves, which guarantees that some favorable properties of the Kleinman-Newton iterates are automatically preserved. Our paper [10] extends line search algorithms and their analyses to inexact solves. An inexact Kleinman-Newton method without line search is analyzed in [13], but some assumptions made in [13] do not hold when low-rank methods are applied to solve the Lyapunov equation iteratively. We extended the inexact Kleinman-Newton method and analysis to integrate the efficient low-rank ADI solver in [10]. This section reviews the main algorithmic components of [10] applied to the projected GCARE (2.15). The following Section 4 then carefully exploits the special structure of the projected GCARE (2.15) for an efficient numerical realization.
At its core our method is an inexact Newton method applied to the GCARE R(X ) = 0, where R(X ) is the Riccati residual (2.16). Given an approximate solution X (k) ∈ R (n v -n p )×(n v -n p ) and a so-called forcing parameter η k ∈ (0, 1), we compute a
Then we compute a step size ξ k ∈ (0, 1] such that the sufficient decrease condition
is satisfied, where β > 0 is a given parameter. The new iterate is
(3.3)
We will discuss below how we compute an S (k) that satisfies (3.1). As we have shown in [10], if the forcing parameters in (3.1) are limited by
then the sufficient decrease condition (3.2) is satisfied for all step sizes ξ k
To ensure convergence of the sequence of iterates {X (k) }, the step sizes ξ k also need to be bounded away from zero. We will state the precise convergence result later, see Theorem 6 below. We use the Armijo rule to compute the step sizes ξ k . This step size rule and others are discussed in [10], as well as conditions that ensure ξ k ≥ ξ min > 0 for all k. Instead of computing the new iterate S (k) as an approximate solution of R (X (k) )S (k) = -R(X (k) ), it is more favorable for our purposes to compute
as an approximate solution of R (X (k) ) (k) are Lyapunov equations, but the right hand side of the latter equation,
where K (k) is defined in (2.18), is low-rank and this will allow the application of the efficient low-rank ADI method (discussed in the next section) to compute X (k+1) . Note that
where
We define the projected Lyapunov residual at any X (k+1) by
(3.7)
The inexactness condition (3.1) means that we have to compute X (k+1) with
such that the corresponding projected Lyapunov residual satisfies
(3.9)
Using the definition (2.16), (3.5), and (3.8), the residual of the projected CARE at (3.3) can be written as
which can be evaluated efficiently for any ξ k , and therefore can be used to efficiently compute a step size ξ k > 0 that satisfies (3.2). The inexact Kleinman-Newton method with line search is summarized in Algorithm 1 below.
Algorithm 1 Inexact Kleinman-Newton method with line search.
Input: A, M, B, C, tol Newton , initial stabilizing iterate X (0) , η ∈ (0, 1), and β ∈ (0, 1 -η) Output: Approximate unique stabilizing solution X ( * ) of GCARE (2.15)
Select η k ∈ (0, η].
Compute X (k+1) that solves the inexact Lyapunov equation
The following convergence theorem for the iterates generated by Algorithm 1 is adopted from [10,Thm. 10] to match the notation of the projected Riccati equation (2.15). Theorem 6. Let (A, B; M) be stabilizable, let (C, A; M) be detectable and assume that for all k, there exists a symmetric positive semi-definite X (k+1) such that (3.8) and (3.9) hold. Furthermore, let X (k) be the iterates generated by Algorithm 1 andA
(i) If the step sizes are bounded away from zero, i.e., ξ k ≥ ξ min > 0 for all k, then R(X (k) ) F → 0. (ii) If in addition the pencils (A (k) , M) are stable for k ≥ k 0 , and X (k) 0 for all k ≥ k 0 , then X (k) → X ( * ) , where X ( * ) 0 is the unique stabilizing solution of the GCARE (2.15).
The main expense in the inexact Kleinman-Newton Algorithm 1 is in Step 3.11b. We apply the real low-rank ADI method, which is detailed in [10] and in [27,Sec. 6.3.1]. This method generates a low-rank approximate solution X (k+1) of the Lyapunov equation in factored form. Compared to the original version of the ADI method [7,20], which is also the version used in Bänsch et al. [4], we use two important modifications of the original ADI method. The first reorganizes the computation to obtain a low-rank representation of the Lyapunov residual in the ADI iterations [8], and the second exploits the fact that the ADI shifts must occur either as a real number or as a pair of complex conjugate numbers to write almost allfoot_0 matrices in the ADI iterations as real matrices [8]. Most importantly, the improved method generates real matrices Z and W , each with few columns, such that ZZ T = X (k+1) satisfies (3.11a) and the corresponding Lyapunov residual L( X (k+1) ) = W W T obeys (3.11). We refer to [10] or [27,Sec. 6.3.1] for details on the derivation of the real low-rank ADI method. The detailed listing of this method is given in Algorithm 2 below.
Algorithm 2 Generalized real-valued low-rank residual ADI method.
W -1
if Im (q ) = 0 then 5:
√ -Re (q ), δ = Re (q ) / Im (q )
9:
10:
= + 1 15: end while Algorithms 1 and 2 work with the projected matrices, but need to be implemented operating on the matrices M,A,B,C. This transformation will be described in the next section.
The inexact Kleinman-Newton Algorithm 1 and the improved ADI Algorithm 2 are derived and stated in terms of the projected matrices in (2.12). As stated before, these matrices are dense, expensive to compute with and the explicit use of the projection needs to be avoided. As before, we use calligraphic font, like X (k) , to denote projected quantities, and roman font, like X (k) , to denote the corresponding quantities without projection.
Regarding the transformation of the iterates in the inexact Kleinman-Newton Algorithm 1, we already know from (2.19) that
To undo the projections, we multiply the Lyapunov equations and the Riccati residuals from the left by l and from the right by T l and replace Steps 6 and 7 in Algorithm 1 by the following.
6: Compute X (k+1) that solves the inexact Lyapunov equation The reason for multiplying by l and T l is that the projection Π emerges. In fact, using (2.12), (2.9), and (4.1), the left hand side in (4.2a) becomes
where
Although the projection Π emerges in (4.3), it will not be computed and used explicitly. We outline the main ideas in the following subsections.
Recall (2.9) and (2.12). We have
To transform the matrices in the improved ADI Algorithm 2 we set
Using (2.12) and (4.1), the linear system in Step 3 of Algorithm 2 is transformed into
We define
(4.7a)
From (2.9) it follows that
(4.7b)
Finally, multiplying (4.6) by l from the left, using (2.9), (4.7a) and (4.7b), the linear system in Step 3 of Algorithm 2 is written as
As it is shown by Heinkenschloss et al. [16] and Bänsch et al. [4] the solution of the projected system (4.8) is equivalent to the solution of the 2 × 2 block system
where " * " indicates that the second block of the solution matrix is not needed. Finally, since K (k) B T is dense, the matrix in (4.9) is written as a low-rank perturbation
and the solution of (4.9) is computed using the Sherman-Morrison-Woodbury formula. See Bänsch et al. [4] or Weichelt [27, p. 67]. We use (2.9), (4.5) to write the projected Lyapunov residual
(4.10)
Rather than computing W and then multiplying by Π , we can update W = Π W directly. In fact, multiplying line 5 in Algorithm 2 with l from the left and using (4.5), (4.7a) yields
where in the last step we have used the M-orthogonality of Π , i.e., Π M = MΠ T and (4.7b). Thus, the projected low-rank residual factor can be accumulated via
without using any explicit projections. Only the initial right hand side W (k) needs to be projected to define
(4.12)
This one projection at the beginning of the ADI method is computed by first solving
(again, " * " indicates that the second block of the solution is not used) and then setting
See Heinkenschloss et al. [16] or Weichelt [27,Lemma 4.1]. This projection is less expensive than a single ADI step and does not considerably increase the overall computation costs. Moreover, the right-hand side W -1 in (4.8), (4.9) can be replaced by W -1 , since
To incorporate this improved ADI method into Algorithm 1, some remaining issues, such as the storage of the Newton step and the projected Riccati residual, need to be addressed. This is done in the next subsection, expanding the statements in [10, Sec. 5.2].
The Newton step S (k) = r S (k) T r is only used in the computation of the step size ξ k , since the inexact Kleinman-Newton step (3.8) directly iterates over the preliminary solution (3.5), and the definition of the feedback matrix in (2.14), this product can be written as
which characterizes the feedback change corresponding to the preliminary or definite new iterate X (k+1) or X (k+1) . Using (2.12), (4.1), and
which characterizes the feedback change corresponding to the preliminary new iterate or X (k+1) or new iterate X (k+1) . Hence, the dense Newton step S (k) is never formed explicitly.
The definition X (k+1) = ZZ T and update in Step 13 of Algorithm 2 implies the formula (k+1) .
The Riccati residual can be written in low-rank form as
This representation can be used to efficiently compute l R(X (k) ) T l F .
In the initial iteration k = 0 with X (0) = 0, (4.18) holds with W (0) = C T and K (k) = 0. Equation (3.10) and L( X
which is of the form (4.18) with
Using (4.5), (4.19) the projected Riccati residual R(X (k+1) ) := l R(X (k+1) ) T l can be written as k+1) . In the second to last equation in (4.20) we have used the identity
which follows from the M-orthogonality of Π and Π T (X (k+1) -
Algorithm 3 Inexact low-rank Kleinman-Newton-ADI for index-2 DAE systems.
Input: M, A, G, B, C , initial feedback K (0) , tol Newton , η ∈ (0, 1), and β ∈ (0, 1 -η)
Compute ADI shifts {q i } nADI i=1 = {q i } nADI i=1 ⊂ C -ordered such that complex pairs form consecutive entries and choose η k ∈ (0, η].
Set
if Im (q ) = 0 then 10:
13: else 14:
15:
16:
= + 1 18:
end while 23:
Compute ξ k ∈ (0, 1) using, e.g., the Armijo rule.
else 26:
end if 28:
31:
where K (0) is an initial stabilizing feedback. Equation (4.20) shows that the Riccati residual R(X (k+1) ) can be computed without any additional explicit projection. The representation (3.10) shows that R(X (k) + ξ k S (k) ) is a quartic polynomial with scalar coefficients. Just as in [10,Sec. 5] this is used for an efficient implementation of the line search computation.
The final feedback at the end of the k + 1-st Newton step is defined via
Only the feedback matrix is needed, but if desired the Riccati iterate can be computed in low-rank form as follows. Assuming the previous Riccati iterate is defined via X (k) = Z (k) Z (k) T and the preliminary solution is defined via X (k+1) = Z (k+1) Z (k+1) T , the new Riccati iterate can be written as
whose size depends on the number of ADI steps in the k-th and k+1-st Newton iteration. The entire process of the inexact low-rank KN-ADI method is depicted in Algorithm 3.
We illustrate the benefits of Algorithm 3 to solve the GCARE associated with the solution of LQR problem (2.1), (2.2) governed by the linearized Navier-Stokes equation. Since our problem set-up is identical to that in the paper by Bänsch et al. [4] and in Weichelt's PhD Thesis [27], we only sketch it here and refer to [4,27] for details. Additional numerical results can be found in [27]. The domain on which the Navier-Stokes and linearized Navier-Stokes equations are posed is shown in Fig. 1. Inflow boundary conditions are posed on the left boundary, no-slip conditions are posed on part of the cylinder boundary and on the top and bottom boundary, and outflow conditions are imposed on the right boundary. Controls are applied on two segments on the cylinder wall (indicated by feed,1 , feed,2 ). Specifically, for each segment a spatial profile is specified, so that the number of inputs in (2.2) is n u = 2. As described in detail in [4,Sec. 2.7], [27,Sec. 4.1.3], an operator is constructed that converts these Dirichlet boundary controls to distributed controls, such that
2). The observations are chosen to be the vertical velocities of the linearized Navier-Stokes equations at the seven points indicated by P obs,i . Thus, y(t) ∈ R 7 , n y = 7. Moreover, we penalize the output by α > 0, i.e., the output equation (2.2c) takes the concrete form y(t) = αCv(t) with C ∈ R 7×n v specified in [4, p. A855], [27,Sec. 4.4.1].
The solution to the steady state Navier-Stokes equation around which is linearized, as well as the linearized Navier-Stokes equations, i.e., the matrices in (2.1), (2.2) are computed using the finite element flow solver NAVIER [1], which uses P 2 -P 1
Taylor-Hood elements and is written in FORTRAN90. The matrices in (2.1), (2.2) are generated using NAVIER and then stored using the so-called matrix market format [12]. The computations for the resulting matrix equations are performed with MATLAB R2012b on a 64-bit CentOS 5.5 server with Intel Xeon X5650 at 2.67 GHz, with 2 CPUs, 12 cores (6 cores per CPU), and 48 GB main memory available.
We conduct experiments with Reynolds number Re = 100, 200, 300, 400, 500, and we use six finite element discretization levels, with Level 1 being the coarsest (shown in Fig. 1). The matrix sizes corresponding to these discretizations are listed in Table 1.
For larger Reynolds number, the matrix pencil (A, M) is not stable (see [4, Fig. 2], [27,Sec. 4.2.3]) and a nonzero initial feedback is needed. We construct the initial feedback K (0) as specified in [4, Sec. 2.7], [27,Sec. 4.2.3]. First, we illustrate the impact of the line search. Fig. 2 shows the convergence of the 'exact' Kleinman-Newton method (i.e., the Lyapunov equation is solved with fixed high residual tolerance) and the inexact Kleinman-Newton method (Algorithm 1 with η k = min{0.1, 0.9 • R(X (k) ) F }) both with and without line search for the LQR problem governed by the discretized linearized Navier-Stokes equations with Re = 500, output weight α = 10 4 , and discretization level 1. There is little difference in the Riccati residuals between the exact and the inexact Kleinman-Newton method. However, there is a big difference between the method with and without line search. Without line search the relative residual grows dramatically in the initial (k = 0) iteration. With line search, the line search is active ξ k < 1 for iterations k = 0, 1, 2 (exact Kleinman-Newton) and iterations k = 0, 1 (inexact Kleinman-Newton). Fig. 2 also shows the Riccati residuals corresponding to ξ k = 1 for the iterations where the line search is active. That the line search is typically only active in the first few iterations has also been observed in other applications of Riccati equations (see, e.g., [5]).
Next, we illustrate the influence of the various improvements to the overall performance of the Algorithm 3. Specifically we compare five set-ups, where 'Setup i' corresponds to a basic version of the Kleinman-Newton-ADI method, and 'Setup v' Finally, adding the line-search in Setup v improves the method further. The number of ADI iterations and linear system solves is reduced by a factor of two. The reduction in ADI iterations also reduced the time for the shift computation. The line search is less than one only in the first iteration and the cost of step size computation is negligible. Comparing the total computation times shows that the algorithm specified in Setup v is 90-times faster than the algorithm specified in Setup i.
As we will see next, the solution of the Riccati equation becomes more difficult as the output weighting α increases. In those cases the speedup of Setup v over Setup i is even more important.
The following numerical tests focus on Algorithm 3, with Setup v. As mentioned before, Algorithm 3 with Setup v will be referred to as 'iKNqLS'. Table 3 documents the performance of iKNqLS applied to our test problem as Reynolds number Re, output weight α, and discretization level change. Table 3a shows that the number of Kleinman-Newton iterations increases moderately with an increasing α and increasing Reynolds number. Similarly, the number of total ADI iterations needed to approximately solve the Lyapunov equations increases with an increasing α and increasing Reynolds number. Furthermore, line search is only necessary for higher Reynolds numbers and higher output weights.
Table 3b shows that for Re ≤ 300, the number of Newton iterations remains nearly constant as the refinement level is increased. For Re ≤ 300 and refinement level greater than two the number of iterations where the step size is less than one is unusually large. We believe that this effect is a result of the instability of the matrix pencil. Solving the first Newton step inexactly might yield an intermediate solution that is slightly (especially in finite precision arithmetic) not stabilizing. Therefore, the following ADI iteration tends to diverge. Our algorithm detects this behavior by monitoring the Riccati and Lyapunov residual continuously. Although this behavior is not covered by the convergence proof in Theorem 6, where a stabilizing solution for k ≥ k 0 is required, our algorithm handles this situation by deleting the last ADI step and performing a line search. This yields convergence in all examples we considered. The relative Riccati residual seems to stagnate for a couple of iterations and, hence, an increasing amount of line search runs is required.
Convergence theory for the exact Kleinman-Newton method guarantees that the matrix pencils are stable if the initial matrix pencil is stable. Thus, another approach to circumvent the appearance of a possibly unstable pencil arising from inexact Lyapunov equation solution is to use a smaller fixed ADI tolerance for the first Newton step. Rather than using the Lyapunov residual tolerance tol ADI = η k R(X (k+1) ) , we set tol ADI = 10 -2 for the first two Newton iterations. If the relative Riccati residual decreases and drops below 5 •10 -1 , the method switches to the iKNqLS scheme (i.e., tol ADI = η k R(X (k+1) ) ). This is referred to as 'exact' start. Using the tol ADI = η k R(X (k+1) ) in all iterations is referred to as inexact start. Table 4 compares both starting procedures for Re ≥ 300 and refinements Level 3-6. The 'exact' start prevents the stagnation of the relative Riccati residual and reduces the number of Newton iterations. The line search is used in at most one iteration. However, the 'exact' solves in the first Newton iterations increase the number of ADI iterations. Therefore, in most cases a decrease in the number of Newton iterations does not translate into a significant decrease in the total number of ADI iterations (and therefore significant decrease in overall computing time) when the 'exact' start is used.
Overall iKNqLS is able to solve the Riccati equation in all cases. Although there is no theoretical justification, our numerics indicate that the inclusion of line search and computationally inexpensive monitoring of the low-rank Riccati and Lyapunov residuals enables the algorithm to successfully cope with intermediate iterates that are nearly not stabilizing.
We have extended our inexact Kleinman-Newton method low-rank ADI solver and line search from [10] to Riccati equations governed by Hessenberg index-2 DAEs. Using the projection idea from Heinkenschloss et al. [16] and Bänsch, Benner [2] we transform the problem governed by the DAE into a 'classical' problem governed by an ODE. Our algorithm in [10] is then applied to this transformed problem. However, the projected ODE is never computed in practice. Instead, a careful exploitation of the problem structure allows the formulation of the algorithm in the original DAE context. We have demonstrated the performance of our Riccati solver to a problem arising in feedback stabilization of Navier-Stokes flow around a cylinder. The numerical results document the impact of various algorithmic components on the overall performance. The algorithmic improvements in this paper lead to approximately 90-times speed-up over a previously used Kleinman-Newton-ADI method. Moreover, we have explored the performance of the new algorithm for various Reynolds numbers, mesh refinement levels and output weights. The new algorithm was able to solve all instances. Moreover, although there is no theoretical justification, our numerics indicate that the inclusion of line search and computationally inexpensive monitoring of the low-rank Riccati and Lyapunov residuals enables the new algorithm to successfully cope with intermediate iterates that are (slightly) not stabilizing.
The linear system solve still has a complex coefficient matrix and thus the intermediate V is complex. This can be avoided along the lines of[18, Remark
4.4], but is not done in our implementation.