1 Introduction

We consider the incompressible Navier–Stokes equations with inputs and outputs

$$\begin{aligned} \dot{v}&= -(v\cdot \nabla )v + \frac{1}{\mathsf {Re}} \varDelta v - \nabla p + \mathcal {B}u, \end{aligned}$$
(1a)
$$\begin{aligned} 0&= \nabla \cdot v, \end{aligned}$$
(1b)
$$\begin{aligned} y&= \mathcal {C}v, \end{aligned}$$
(1c)

and the question when a linear output-feedback controller \(\mathcal {K}:y \mapsto u\) can stabilize this nonlinear system around a possibly unstable steady state in the presence of system uncertainties. Here, the variables v and p describe the evolution of the velocity and pressure fields in a given flow setup that is parametrized through the Reynolds number \(\mathsf {Re}\), the operator \(\mathcal {B}\) models the actuation through the controls, and \(\mathcal {C}\) is the output operator.

We will approach this question through a semi-discrete and linearized approximation to (1), model order reduction to cope with the high dimensionality of the controller design problem, and the design of controllers that can compensate for a large class of system uncertainties. Basically, our argument is that discretization and model reduction errors are of the same nature such that a proven robustness margin can possibly overcome unmodeled uncertainties, too. Anyways, in order to potentially work in physical setups, any model-based controller needs a certain robustness against inevitable model errors. This rules out the standard linear quadratic Gaussian (LQG) design that has no guaranteed stability margin [21]. A general remedy is provided by \(\mathcal {H}_{\infty }\)-controllers that, provably, can compensate for linearization errors [7, 14, 30], discretization errors [8, 19], and truncation errors [37].

The \(\mathcal {H}_{\infty }\)-theory roots in the 1980s [48]; see also [23, 24] for the historical background. In view of its application in simulations, the development of state-space formulations [22, 24] meant a breakthrough since it came with general formulas for the controller design based on the solutions of indefinite Riccati equations. Nonetheless, the computational effort for solving these Riccati equations is significant so that, up to now, this design approach has rarely been considered in large-scale simulations let alone the case where algebraic constraints are present, that is for differential-algebraic equations (DAEs) or descriptor systems. If one leaves aside the offline effort for designing the controller, the theory seems well suited for large-scale problems since the controllers allow for a low-dimensional approximation with a-priori estimates on performance and robustness [37]. The involved reduction is based on balanced truncation and the related \({\mathrm {LQG}}\)-approach has been investigated for descriptor systems in [36]. Thus, with flow control in mind, the solution of large-scale Riccati equations related to DAEs enables the use of the general \(\mathcal {H}_{\infty }\)-theory. It has been acknowledged that for general DAE systems, the standard symmetric Riccati equations are only suitable under very restrictive conditions [6]. A nonsymmetric version has been shown to provide a true generalization for the \(\mathcal {H}_{\infty }\)-controller [10] to the index-1 case and is widely applicable for LQG-design for impulse controllable descriptor systems [31, 45].

In this work, we consider the incompressible Navier–Stokes equations with control inputs in the momentum equation and velocity measurements only, so that the input-to-output behavior can be equivalently realized as a system of ordinary differential equations (ODEs). Still, we keep the DAE structure since the corresponding transformation will not be available in practice. Then, the challenge is to realize what is suggested by the ODE theory without explicitly resorting to transformations and projections. Therefor, we adapt the established technique of realizing the projections through the solution of saddle-point systems [2, 28, 32, 46], such that structure and sparsity are preserved during all operations.

The \(\mathcal {H}_{\infty }\)-feedback control for the 2D incompressible Navier–Stokes equations has been treated in [20] from a theoretical perspective. Therein, unmodeled boundary inputs are considered and existence of optimal feedback solutions is shown with the help of Riccati equations. The general \(\mathcal {H}_{\infty }\)-control problem is discussed in [3, Ch. 5]. Here, summing up the relevant research, the established state-space theory is adapted to the infinite-dimensional incompressible Navier–Stokes equations.

In this paper, we explore the Riccati-based \(\mathcal {H}_{\infty }\)-controller design for spatially discretized two-dimensional incompressible flows. To this end,

  • we adapt the theory for the implicit treatment of the incompressibility constraint to the \(\mathcal {H}_{\infty }\)-optimization problem,

  • we leverage \(\mathcal {H}_{\infty }\)-balanced truncation to reduce the dimension of the controller design problem,

  • we provide numerically accessible formulas for a-priori estimation of the robustness of the controller with respect to the \(\mathcal {H}_{\infty }\)-balanced truncation model reduction error as well as linearization errors, and

  • illustrate the performance in two challenging numerical examples.

As a result, we provide a complete numerical approach that makes \(\mathcal {H}_{\infty }\)-controller design feasible for large-scale Navier–Stokes systems and provides computable bounds on the robustness of the performance with respect to both linearization errors [14] and model reduction errors [37]. The situation that the controller is based on inexact linearizations is relevant in applications and fits well into the presented framework; cf. [8, 14, 30]. The presented numerical studies are based on the popular setup of the two-dimensional wake behind a cylinder; see, for examples, [17, 26, 27, 39, 47].

This paper is organized as follows: In Sect. 2, we introduce the concepts of \(\mathcal {H}_{\infty }\)-controller design and truncation via the solution of Riccati equations and how the theory extends to semi-discretized incompressible Navier–Stokes equations. In Sect. 3, we discuss numerical methods for the solution of large-scale Riccati equations that implicitly respect the incompressibility constraint and provide a summary of steps for the design of robust low-dimensional controllers with the accompanying formulas. By means of two flow setups, we report on the performance of the resulting low-dimensional controllers in Sect. 4. The paper is concluded in Sect. 5.

2 Mathematical basics

In this section, we introduce the basic concept and state-space approach of robust \(\mathcal {H}_{\infty }\)-controller design and how the relevant formulas are realized for incompressible flow control setups.

2.1 Riccati-based \(\mathcal {H}_{\infty }\) controller design

For practical applications, the design of controllers that provide a certain robustness against disturbances is essential, as neither models nor numerical computations are fully able to match reality. The \(\mathcal {H}_{\infty }\)-control theory provides the design of such controllers. In general, linear time-invariant systems of the form

$$\begin{aligned} E {\dot{x}}(t)&= A x(t) + B_{1} w(t) + B_{2} u(t), \end{aligned}$$
(2a)
$$\begin{aligned} z(t)&= C_{1} x(t) + D_{11} w(t) + D_{12} u(t), \end{aligned}$$
(2b)
$$\begin{aligned} y(t)&= C_{2} x(t) + D_{21} w(t) + D_{22} u(t), \end{aligned}$$
(2c)

with \(E, A \in \mathbb {R}^{n \times n}\), \(B_{1} \in \mathbb {R}^{n \times m_{1}}\), \(B_{2} \in \mathbb {R}^{n \times m_{2}}\), \(C_{1} \in \mathbb {R}^{p_{1} \times n}\), \(C_{2} \in \mathbb {R}^{p_{2} \times n}\), \(D_{11} \in \mathbb {R}^{p_{1} \times m_{1}}\), \(D_{12} \in \mathbb {R}^{p_{1} \times m_{2}}\), \(D_{21} \in \mathbb {R}^{p_{2} \times m_{1}}\), and \(D_{22} \in \mathbb {R}^{p_{2} \times m_{2}}\), are considered. In (2), the internal states x are influenced by the control inputs u and disturbances w, and the user can observe the measurements y and performance outputs z of the system. The basic aim is to find a feedback controller \(K:y \mapsto u\) that internally stabilizes (2). Considering (2) in frequency domain allows the formulation of the system’s input-to-output relation in terms of its (partitioned with respect to the different inputs and outputs) transfer function

$$\begin{aligned} G(s) = \left[ \begin{array}{l} C_{1} \\ C_{2} \end{array}\right] (sE - A)^{-1} \left[ \begin{array}{ll} B_{1}&B_{2} \end{array}\right] + \left[ \begin{array}{ll} D_{11} &{} D_{12} \\ D_{21} &{} D_{22} \end{array}\right] =: \left[ \begin{array}{ll} G_{11}(s) &{} G_{12}(s) \\ G_{21}(s) &{} G_{22}(s) \end{array}\right] . \end{aligned}$$

Let also K(s) denote the transfer function of the controller K, the disturbance-to-performance behavior of the system can be formulated as

$$\begin{aligned} \begin{aligned} Z(s)&= \Big (G_{11}(s) + G_{12}(s) K(s) \big (I_{p_{2}} - G_{22}(s) K(s)\big )^{-1} G_{21}(s) \Big ) W(s) \\&=: \mathcal {F}(G, K) W(s), \end{aligned} \end{aligned}$$
(3)

where Z and W are the Laplace transforms of the performances and disturbances, and \(I_{p_{2}}\) denotes the \(p_{2}\)-dimensional identity matrix. \(\mathcal {F}\) is a lower linear fractional transformation. With (3), the optimal \(\mathcal {H}_{\infty }\)-control problem is to find a controller K that solves

$$\begin{aligned} \min \limits _{K\,\text {stabilizing}} \Vert \mathcal {F}(G, K) \Vert _{\mathcal {H}_{\infty }} =: \gamma _{\mathsf {opt}}, \end{aligned}$$
(4)

where \(\Vert .\Vert _{\mathcal {H}_{\infty }}\) is the \(\mathcal {H}_{\infty }\)-norm. In general, this optimization problem (4) is too difficult to solve. Instead, one can consider a relaxation in terms of the suboptimal \(\mathcal {H}_{\infty }\)-control problem: Find a stabilizing controller K such that

$$\begin{aligned} \Vert \mathcal {F}(G, K) \Vert _{\mathcal {H}_{\infty }} < \gamma , \end{aligned}$$
(5)

which is then solved successively for decreasing robustness margins \(\gamma \rightarrow \gamma _{\mathsf {opt}}\).

A state-space theory has been developed that provides formulas for such suboptimal controllers. For ease of notation and theoretical derivations, we use the following set of assumptions:

Assumption 1

(Simplifying and necessary assumptions for the \(\mathcal {H}_{\infty }\) -controller design) In the formulation of system (2), we have that

  1. 1.

    E is invertible,

  2. 2.

    \(D_{11}, D_{22} = 0\),

  3. 3.

    \((sE - A, B_{1})\) is stabilizable and \((sE - A, C_{1})\) is detectable,

  4. 4.

    \((sE - A, B_{2})\) is stabilizable and \((sE - A, C_{2})\) is detectable,

  5. 5.

    \(D_{12}^{\mathsf {T}} \left[ \begin{array}{ll} C_{1}&D_{12} \end{array}\right] = \left[ \begin{array}{ll} 0&I \end{array}\right]\), and \(\left[ \begin{array}{ll} B_{1} \\ D_{21} \end{array}\right] D_{21}^{\mathsf {T}} = \left[ \begin{array}{ll} 0 \\ I \end{array}\right]\).

We will comment on the assumptions in Remark 1 after we have introduced the procedure to design an \(\mathcal {H}_{\infty }\)-robust controller.

Proposition 1

(Existence of \(\mathcal {H}_{\infty }\)-controllers [22, 49]) Given a \(\gamma > \gamma _{\mathsf {opt}}\), there exists an admissible controller K if and only if there are unique, symmetric positive semi-definite stabilizing solutions \(X_{\mathcal {H}_{\infty }}\) and \(Y_{\mathcal {H}_{\infty }}\) to the regulator and filter \(\mathcal {H}_{\infty }\) -Riccati equations

$$\begin{aligned} A^{\mathsf {T}} X E + E^{\mathsf {T}} X A - E^{\mathsf {T}} X (B_{2}B_{2}^{\mathsf {T}}-\gamma ^{-2}B_{1}B_{1}^{\mathsf {T}} ) X E + C_{1}^{\mathsf {T}} C_{1}&= 0, \end{aligned}$$
(6a)
$$\begin{aligned} A Y E^{\mathsf {T}} + E Y A^{\mathsf {T}} - E Y (C_{2}^{\mathsf {T}} C_{2} - \gamma ^{-2} C_{1}^{\mathsf {T}} C_{1} )YE^{\mathsf {T}} + B_{1} B_{1}^{\mathsf {T}}&= 0, \end{aligned}$$
(6b)

and, additionally,

$$\begin{aligned} \gamma ^{2} > \lambda _{\max }(Y_{\mathcal {H}_{\infty }} E^{\mathsf {T}} X_{\mathcal {H}_{\infty }} E), \end{aligned}$$
(7)

where \(\lambda _{\max }(M)\) denotes the maximum eigenvalue of the matrix M.

In the case that Proposition 1 holds, a stabilizing controller solving (5), known in the literature as the central or minimum entropy controller, is given via

$$\begin{aligned} K:\left\{ \begin{array}{l} \widetilde{E}{\dot{\tilde{x}}}(t) = \widetilde{A}\tilde{x}(t) + \widetilde{B}y(t),\\ u(t) = \widetilde{C}\tilde{x}(t), \end{array}\right. \end{aligned}$$

where the system matrices can be computed as

$$\begin{aligned} \widetilde{E}&= E, \end{aligned}$$
(8a)
$$\begin{aligned} \widetilde{A}&= A + E Y_{\mathcal {H}_{\infty }}(\gamma ^{-2} C_{1}^{\mathsf {T}}C_{1} - C_{2}^{\mathsf {T}}C_{2}) - B_{2}B_{2}^{\mathsf {T}} X_{\mathcal {H}_{\infty }} E Z_{\mathcal {H}_{\infty }}, \end{aligned}$$
(8b)
$$\begin{aligned} \widetilde{B}&= E Y_{\mathcal {H}_{\infty }} C_{2}^{\mathsf {T}}, \end{aligned}$$
(8c)
$$\begin{aligned} \widetilde{C}&= -B_{2}^{\mathsf {T}} X_{\mathcal {H}_{\infty }} E Z_{\mathcal {H}_{\infty }}, \end{aligned}$$
(8d)

with \(Z_{\mathcal {H}_{\infty }} = (I_{n} - \gamma ^{-2}Y_{\mathcal {H}_{\infty }} E^{\mathsf {T}} X_{\mathcal {H}_{\infty }} E)^{-1}\); see [22, Sec. III.C] for the general case and also [37, Eqs. (16) and (17)] for the normalized case described in the following section.

Remark 1

(Comments on the simplifying and necessary assumptions) With E invertible, our Assumption 1 is equivalent to the standard set of assumptions made in \(\mathcal {H}_{\infty }\)-robust controller design; cf., e.g., [22, Sec. II.A].

The assumptions on stabilizability and detectability with respect to \(B_2\) and \(C_2\) are clearly necessary for the existence of a stabilizing output feedback controller in the case of no disturbances. Stabilizability and detectability with respect to \(B_1\) and \(C_1\) are assumed to simplify both the formulas and the derivations, and can be relaxed towards conditions that ensure that the Riccati equations (6) are solvable in the \(\mathcal {H}_2\) case, i.e., as \(\gamma \rightarrow \infty\). A further relaxation would possibly rule out the proposed Riccati-based approach; cf. [49, Ch. 17.1].49, Ch. 17.1].

Assumption 1 (2) is basically made for simplicity of the formulas. With the transformation \(K = (I + D_{22}{\hat{K}})^{-1}{\hat{K}}\) that maps a controller for \(D_{22}=0\) onto an equivalent controller for the case that \(D_{22}\ne 0\), the assumption of \(D_{22}\) poses no restriction at all. Similarly, the restriction to \(D_{11}= 0\) can be lifted but complicates the formulas a lot.

Finally, the necessary part of Assumption 1 (5) is that these matrices have full rank. The explicit structure can then be achieved by input and output transformations. If, however, the full-rank condition is not met, then the formulation may lead to a singular control problem; see, e.g., [49 , Ch. 17.1].

2.2 The normalized \(\mathcal {H}_{\infty }\) problem and low-rank robust controllers

The target application of this work is output-based feedback control of a nonlinear system that is robust against system uncertainties stemming from linearization and reduction errors. As a general model for the linearization uncertainty, we will assume that the system matrix A is subjected to an additive perturbation. In this case, in a feedback arrangement, disturbance inputs will be induced by the measurements of the perturbed state and enter the system as a perturbation of the control input. Accordingly, in (2), one can consider \(B_{1} = B_{2} =:B\) and \(C_{1} = C_{2} =: C\), which leads to the so-called normalized \(\mathcal {H}_{\infty }\)-control problem; see [37].

A robust \(\mathcal {H}_{\infty }\)-controller based on the normalized problem then applies in our context as follows. First, its robustness margin \(\gamma\) can be weighed up against the system error that arises from the truncation of certain states of the controller. Since a full-order controller would be of the same size as the system, such a truncation significantly supports the efficient evaluation of the feedback law during a simulation. Second, the robustness can cover the system error that comes from an inaccurate linearization. In this section, we discuss both application scenarios and provide relevant a-priori and a-posteriori estimates.

As for controller truncation, we consider the \(\mathcal {H}_{\infty }\)-balanced truncation approach as it has been reported in [37]. Therefor, we introduce another special case of (2), namely the so-called normalized LQG system

$$\begin{aligned} E {\dot{x}}(t)&= Ax(t) + B w_{1}(t) + B u(t), \end{aligned}$$
(9a)
$$\begin{aligned} z_{1}(t)&= C x(t), \end{aligned}$$
(9b)
$$\begin{aligned} z_{2}(t)&= u(t), \end{aligned}$$
(9c)
$$\begin{aligned} y(t)&= C x(t) + w_{2}(t). \end{aligned}$$
(9d)

In the terms of the general system (2), the particular structural assumptions made for (9) mean that

$$\begin{aligned} B_2 = B, \quad C_2=C, \quad B_1 = \left[ \begin{array}{ll} B&0 \end{array}\right] , \quad C_1 = \left[ \begin{array}{l} C \\ 0 \end{array}\right] \end{aligned}$$

and

$$\begin{aligned} D_{11} = 0, \quad D_{22}=0, \quad D_{21} = \left[ \begin{array}{ll} 0&I \end{array}\right] , \quad D_{12} = \left[ \begin{array}{l} 0 \\ I \end{array}\right] , \end{aligned}$$

and imply that Assumption 1 is fulfilled if \((sE-A, B)\) is stabilizable and \((sE-A, C)\) is detectable. Moreover, this form of \(B_1\), \(B_2\), \(C_1\) and \(C_2\) implies that

$$\begin{aligned} B_{1} B_{1}^{\mathsf {T}}= B_{2} B_{2}^{\mathsf {T}} = B B^{\mathsf {T}} \quad \text {and} \quad C_{1}^{\mathsf {T}} C_{1} = C_{2}^{\mathsf {T}} C_{2} = C^{\mathsf {T}} C \end{aligned}$$

and, from Proposition 1, that the existence of an admissible controller, which solves the suboptimal \(\mathcal {H}_{\infty }\)-control problem (5) for the gain \(\gamma\), is equivalent to the existence of symmetric positive semi-definite matrices \(X_{\mathcal {H}_{\infty }}\) and \(Y_{\mathcal {H}_{\infty }}\), which solve

$$\begin{aligned} A Y_{\mathcal {H}_{\infty }} E^{\mathsf {T}} + E Y_{\mathcal {H}_{\infty }} A^{\mathsf {T}} - (1 - \gamma ^{-2}) E Y_{\mathcal {H}_{\infty }} C ^{\mathsf {T}}C Y_{\mathcal {H}_{\infty }} E^{\mathsf {T}} + B B ^{\mathsf {T}}&= 0, \end{aligned}$$
(10a)
$$\begin{aligned} A^{\mathsf {T}} X_{\mathcal {H}_{\infty }} E + E^{\mathsf {T}} X_{\mathcal {H}_{\infty }} A - (1 - \gamma ^{-2}) E^{\mathsf {T}} X_{\mathcal {H}_{\infty }} B B ^{\mathsf {T}} X_{\mathcal {H}_{\infty }} E + C^{\mathsf {T}} C&= 0, \end{aligned}$$
(10b)

such that the spectrum condition (7) holds and the matrix pencils \(sE - (A - (1 - \gamma ^{-2}) E Y_{\mathcal {H}_{\infty }} C ^{\mathsf {T}}C)\) and \(sE - (A - (1 - \gamma ^{-2}) B B ^{\mathsf {T}} X_{\mathcal {H}_{\infty }} E)\) are stable.

By means of the two matrices \(X_{\mathcal {H}_{\infty }}\) and \(Y_{\mathcal {H}_{\infty }}\), the so-called \(\mathcal {H}_{\infty }\)-balanced truncation of the system can be computed [37, Prop. 4.10] that enables the truncation of system or controller states with an a-priori control on the approximation error. A practical implementation of this \(\mathcal {H}_{\infty }\)-balanced truncation (HINFBT) for large-scale systems is shown in Algorithm 1 that uses the square root balancing approach with approximating low-rank factorizations of \(X_{\mathcal {H}_{\infty }}\) and \(Y_{\mathcal {H}_{\infty }}\). An implementation of Algorithm 1 for the dense system case can be found in [11].

An error bound for the approximation computed by Algorithm 1 is given in terms of normalized coprime factorizations:

Proposition 2

(\(\mathcal {H}_{\infty }\)-balanced truncation error bound [37]) Let \(G(s) = M(s)^{-1} N(s)\) and \(\widehat{G}(s) = \widehat{M}(s)^{-1} \widehat{N}(s)\) be normalized left coprime factorizations of the original system G and the reduced-order model \(\widehat{G}\) obtained by HINFBT (Algorithm 1) with the robustness margin \(\gamma > \gamma _{\mathsf {opt}}\), respectively. Then, a bound for the approximation error is given by

$$\begin{aligned} \left\Vert \left[ \begin{array}{ll} \beta (N - \widehat{N})&\; M - \widehat{M}\end{array}\right] \right\Vert _{\mathcal {H}_{\infty }}&\le 2 \sum \limits _{k = r + 1}^{n} \frac{\beta \sigma _{k}}{\sqrt{1 + \beta ^{2} \sigma _{k}^{2}}}, \end{aligned}$$
(11)

with \(\beta = \sqrt{1 - \gamma ^{-2}}\) and the characteristic \(\mathcal {H}_{\infty }\) -values \(\sigma _{k}\).

figure a

An important consequence of Proposition 2 is that the robustness margin \(\gamma\) can be put in context with the approximation error of the reduced-order model and, consequently, the final reduced-order controller. Thereby, one can estimate the size of the reduced-order controller, which is needed to still stabilize the original system.

In general, let \(G(s) = M(s)^{-1} N(s)\) and \(\widehat{G}(s) = \widehat{M}(s)^{-1} \widehat{N}(s)\) be as in Proposition 2 the normalized left coprime factorizations of the original system and the reduced-order approximation computed by HINFBT (Algorithm 1), and define

$$\begin{aligned} \beta \hat{\epsilon }:= \left\Vert \left[ \begin{array}{ll} \beta (N - \widehat{N})&\; M - \widehat{M}\end{array}\right] \right\Vert _{\mathcal {H}_{\infty }} \end{aligned}$$
(12)

to be the scaled coprime factor error, with \(\beta = \sqrt{1 - \gamma ^{-2}}\) and \(\gamma\) used in HINFBT (Algorithm 1). Then, it has been shown in [37, Cor. 5.5] that a sufficient condition for the reduced-order central controller \(\widehat{K}\) based on \(\widehat{G}\) to stabilize the original system G is

$$\begin{aligned} \hat{\epsilon }(\beta + \hat{\gamma }) < 1, \end{aligned}$$
(13)

where \(\hat{\gamma }:= \Vert \mathcal {F}(\widehat{G}, \widehat{K}) \Vert _{\mathcal {H}_{\infty }}\). The condition given in (13) is not practical since it can only be evaluated after the computation of the reduced-order model and the computation of the \(\mathcal {H}_{\infty }\)-norm of a large-scale error system is needed. A more practical alternative to (13) can be derived from the a-priori estimate in Proposition 2. Consider the HINFBT error bound (11) to define

$$\begin{aligned} \beta \epsilon := 2 \sum \limits _{k = r + 1}^{n} \frac{\beta \sigma _{k}}{\sqrt{1 + \beta ^{2} \sigma _{k}^{2}}}, \end{aligned}$$

or, more precisely,

$$\begin{aligned} \epsilon := 2 \sum \limits _{k = r + 1}^{n} \frac{\sigma _{k}}{\sqrt{1 + \beta ^{2} \sigma _{k}^{2}}}. \end{aligned}$$
(14)

Note that \(\hat{\epsilon }\le \epsilon\) necessarily holds. Since \(\widehat{G}\) is constructed by HINFBT (Algorithm 1) it also holds that \(\hat{\gamma }\le \gamma\) for \(\gamma\) the robustness margin used in Algorithm 1. Therefor, a sufficient a-priori condition for stabilization of the full-order system by the reduced-order central controller is

$$\begin{aligned} \epsilon (\beta + \gamma ) < 1. \end{aligned}$$
(15)

Remark 2

(Performance loss of reduced-order controller [37]) We need to note that while a reduced-order controller \(\widehat{K}\) guarantees the stabilization of the full-order plant using conditions like (13) or (15), it negatively influences the performance gain of the closed-loop system. One can show that if the full-order central controller (8) satisfies \(\Vert \mathcal {F}(G, K) \Vert _{\mathcal {H}_{\infty }} < \gamma\), then for the reduced-order controller \(\widehat{K}\) based on a HINFBT reduced-order model (Algorithm 1) using the performance gain \(\gamma\) it holds that

$$\begin{aligned} \Vert \mathcal {F}(G, \widehat{K}) \Vert _{\mathcal {H}_{\infty }} \le \hat{\gamma }+ \frac{\hat{\epsilon }(1 + \hat{\gamma })(1 + \beta + \hat{\gamma })}{1 - \hat{\epsilon }(\beta + \hat{\gamma })} < \gamma + \frac{\epsilon (1 + \gamma )(1 + \beta + \gamma )}{1 - \epsilon (\beta + \gamma )}, \end{aligned}$$
(16)

where \(\hat{\epsilon }\) and \(\epsilon\) describe the approximation error and its a-priori bound as in (12) and (14), \(\hat{\gamma }= \Vert \mathcal {F}(\widehat{G}, \widehat{K}) \Vert _{\mathcal {H}_{\infty }} \le \gamma\) and \(\beta = \sqrt{1 - \gamma ^{-2}}\). In practice, the bounds in (16) quickly converge to \(\hat{\gamma }\) or \(\gamma\), respectively, for increasing reduced order.

Another purpose of the robustness margin \(\gamma\) is to ensure stability of the closed-loop system in the presence of system uncertainties. The amount of uncertainty that is guaranteed to be compensated is given in the following proposition.

Proposition 3

(Stabilization of disturbed systems [35, Cor. 3.7]) Given the system \(G(s) = M(s)^{-1} N(s)\) of the form (9) and a stabilizing controller K that solves the sub-optimal \(\mathcal {H}_{\infty }\) -controller problem (5), i.e., \(\Vert \mathcal {F}(G, K) \Vert _{\mathcal {H}_{\infty }} < \gamma\) holds for a given robustness margin \(\gamma\). Let \(G_{\varDelta }(s) = M_{\varDelta }(s)^{-1} N_{\varDelta }(s)\) be another system in normalized LQG form (9), then K is guaranteed to also stabilize \(G_{\varDelta }\) if

$$\begin{aligned} \left\Vert \left[ \begin{array}{ll} N - N_{\varDelta }&\; M - M_{\varDelta } \end{array}\right] \right\Vert _{\mathcal {H}_{\infty }} < \gamma ^{-1}. \end{aligned}$$
(17)

holds.

In view of stabilization of incompressible Navier–Stokes equations by linear output-feedback controllers, the following considerations are relevant with respect to the estimate (17). It has been shown that an error in the linearization used for controller design, smoothly transfers to a coprime factor perturbation in the transfer function; see [7, 30]. Accordingly, with increasing accuracy in the computation of the linearization, the difference becomes arbitrarily small such that, eventually, a robust controller based on a numerically computed linearization will be able to stabilize the system.

We end this section with a summary of the a-priori-type conditions for the construction of reduced-order stabilizing controllers in form of the following theorem.

Theorem 1

(Sufficient a-priori conditions for stabilization of disturbed systems) Given a system \(G = M(s)^{-1} N(s)\) in normalized LQG form (9) and a reduced-order system \(\widehat{G}\) computed by HINFBT (Algorithm 1) with the robustness margin \(\gamma\). The central controller \(\widehat{K}\) based on \(\widehat{G}\) is guaranteed to stabilize the full-order system G if

$$\begin{aligned} \epsilon (\beta + \gamma ) < 1, \end{aligned}$$

where \(\beta = \sqrt{1 - \gamma ^{-2}}\), and \(\epsilon\) is computed from the truncated characteristic \(\mathcal {H}_{\infty }\) -values (14). Additionally, the reduced-order central controller \(\widehat{K}\) is guaranteed to stabilize all systems \(G_{\varDelta }(s) = M_{\varDelta }(s)^{-1} N_{\varDelta }(s)\) for which

$$\begin{aligned} \left\Vert \left[ \begin{array}{ll} N - N_{\varDelta }&\; M - M_{\varDelta } \end{array}\right] \right\Vert _{\mathcal {H}_{\infty }} < \big (\gamma _{\mathrm {G}\widehat{\mathrm {K}}} \big )^{-1} \end{aligned}$$

holds, and where

$$\begin{aligned} \gamma _{\mathrm {G}\widehat{\mathrm {K}}} = \gamma + \frac{\epsilon (1 + \gamma )(1 + \beta +\gamma )}{1 - \epsilon (\beta + \gamma )}. \end{aligned}$$

Note that this theorem allows to adaptively choose the size of the reduced-order controller depending on the stabilization of the original full-order system and the amount of disturbances that shall be compensated by the reduced-order controller.

2.3 Semi-discretization and linearization of Navier–Stokes equations

Next, we briefly discuss how a linear finite dimensional ODE system can be derived as base of the controller design for the nonlinear incompressible Navier–Stokes equations. Details of the discretization and the modeling of boundary control will be discussed together with the numerical examples in Sect. 4.

A finite elements discretization of the incompressible Navier–Stokes equations (1) leads to the semi-discrete system:

$$\begin{aligned} E \dot{v}&= A_{\mathsf {S}}v + F(v) + J^{\mathsf {T}}p + Bu + f, \end{aligned}$$
(18a)
$$\begin{aligned} 0&= Jv + g, \end{aligned}$$
(18b)
$$\begin{aligned} y&= Cv, \end{aligned}$$
(18c)

where \(E \in \mathbb {R}^{n_{\mathrm {v}}\times n_{\mathrm {v}}}\) is the (invertible) mass matrix, \(A_{\mathsf {S}}\in \mathbb {R}^{n_{\mathrm {v}}\times n_{\mathrm {v}}}\) is the discrete approximation of \(\frac{1}{\mathsf {Re}}\varDelta\), \(F:\mathbb {R}^{n_{\mathrm {v}}} \rightarrow \mathbb {R}^{n_{\mathrm {v}}}\) models the convection, \(J \in \mathbb {R}^{n_{\mathrm {p}}\times n_{\mathrm {v}}}\) and \(J^{\mathsf {T}}\) represent the discrete divergence and gradient, \(B \in \mathbb {R}^{n_{\mathrm {v}}\times m}\) and \(C \in \mathbb {R}^{p \times n_{\mathrm {v}}}\) are the discretized input and output operators, and where \(f \in \mathbb {R}^{n_{\mathrm {v}}}\) and \(g \in \mathbb {R}^{n_{\mathrm {p}}}\) are the inhomogeneities that arise from the inclusion of the inflow boundary condition in strong form. Let \(v_\infty\) be the steady state with \(u = 0\) so that with \(v_\delta = v-v_\infty\) and \(A^{(\infty )}:= A_{\mathsf {S}}+ (\partial _v F)(v_\infty )\),

$$\begin{aligned} E {{\dot{v_\delta }}}&= A^{(\infty )}v_\delta + J^{\mathsf {T}}p + Bu, \end{aligned}$$
(19a)
$$\begin{aligned} 0&= Jv_\delta , \end{aligned}$$
(19b)
$$\begin{aligned} y_\delta&= Cv_\delta , \end{aligned}$$
(19c)

provides a linearization that can be used for regulating the deviation \(v_\delta\) from the steady state and, thus, for designing a controller for stabilizing \(v_\infty\). With the assumption that the chosen finite elements scheme is LBB-stable, one can define a discrete realization of the Leray projector

$$\begin{aligned} \varPi ^{\mathsf {T}} := I_{n_{\mathrm {v}}} - E^{-1}J^{\mathsf {T}}(JE^{-1}J^{\mathsf {T}})^{-1}J, \end{aligned}$$
(20)

that maps v(t) into the kernel of J along the orthogonal complement (in the inner product induced by the mass matrix E) of \(J^{\mathsf {T}}\). Making use of \(\varPi\) and the identities \(\varPi E = E \varPi ^{\mathsf {T}}\)—which holds for symmetric E—and \(\varPi ^{\mathsf {T}}v_\delta = v_\delta\), we can eliminate the discrete pressure p and the algebraic constraint \(0 = Jv_\delta\) and rewrite (19) as ODE system

$$\begin{aligned} E {{\dot{v_\delta }}}&= \varPi A^{(\infty )}\varPi ^{\mathsf {T}}v_\delta + \varPi Bu, \end{aligned}$$
(21a)
$$\begin{aligned} y_\delta&= C\varPi ^{\mathsf {T}}v_\delta . \end{aligned}$$
(21b)

For such a realization in terms of an ODE system (21), \(\mathcal {H}_{\infty }\)-robust controllers can be defined via the solutions to the projected \(\mathcal {H}_{\infty }\)-Riccati equations

$$\begin{aligned}&\varPi \big (A^{(\infty )}\big )^{\mathsf {T}}\varPi ^{\mathsf {T}}XE + E^{\mathsf {T}}X \varPi A^{(\infty )}\varPi ^{\mathsf {T}} \nonumber \\&\quad -(1 - \gamma ^{-2}) E^{\mathsf {T}} X \varPi B B^{\mathsf {T}} \varPi ^{\mathsf {T}} X E + \varPi C^{\mathsf {T}} C\varPi ^{\mathsf {T}} = 0, \end{aligned}$$
(22a)
$$\begin{aligned}&\varPi A^{(\infty )}\varPi ^{\mathsf {T}} Y E^{\mathsf {T}} + E Y \varPi \big (A^{(\infty )}\big )^{\mathsf {T}} \varPi ^{\mathsf {T}} \nonumber \\&\quad -(1 - \gamma ^{-2}) E Y\varPi C^{\mathsf {T}}C\varPi ^{\mathsf {T}}YE^{\mathsf {T}} + \varPi BB^{\mathsf {T}}\varPi ^{\mathsf {T}} = 0. \end{aligned}$$
(22b)

In the context of nonlinear flow control, the problem of linearization errors and compensation by reduced-order controllers has been considered in [14]. There it is mentioned that linearization errors in (19) from (18) result in additive disturbances on the transfer function, which can be handled as disturbances in the coprime factorization. In fact, the coprime factorization can be explicitly written down.

Proposition 4

(Normalized left coprime factorization for flow problems [14]) Given the system (19) such that (21) restricted to the image of \(\varPi\) is detectable. The normalized left coprime factorization of (19) is given by

$$\begin{aligned} \left[ \begin{array}{ll} N(s)&M(s) \end{array}\right]&= \mathcal {C}(s \mathcal {E}- \mathcal {A})^{-1} \left[ \begin{array}{ll} \mathcal {B}&-\mathcal {L}\end{array}\right] + \left[ \begin{array}{ll} 0&I_{p} \end{array}\right] , \end{aligned}$$

where the system matrices are

$$\begin{aligned} \begin{aligned} \mathcal {E}&= \left[ \begin{array}{ll} E &{} 0 \\ 0 &{} 0 \end{array}\right] ,&\mathcal {A}&= \left[ \begin{array}{ll} A^{(\infty )}- (1 - \gamma ^{-2}) E Y_{\mathcal {H}_{\infty }}C^{\mathsf {T}}C &{} J^{\mathsf {T}} \\ J &{} 0\end{array}\right] ,\\ \mathcal {C}&= \left[ \begin{array}{ll} C&0 \end{array}\right] ,&\mathcal {B}&= \left[ \begin{array}{ll} B \\ 0 \end{array}\right] ,\quad \mathcal {L}= \left[ \begin{array}{ll} (1 - \gamma ^{-2}) E Y_{\mathcal {H}_{\infty }} C^{\mathsf {T}} \\ 0 \end{array}\right] , \end{aligned} \end{aligned}$$

and with \(Y_{\mathcal {H}_{\infty }}\), the stabilizing solution of (22).

The representation in Proposition 4 allows for the explicit computation of disturbances on the coprime factorizations, e.g., for the evaluation of the necessary robustness as in (17), accessible to numerical simulations.

3 Numerical methods

3.1 Projector-free realization for incompressible flows

In principle, the projected Riccati equations (22) could be treated by established Riccati equation solvers. However, in practice, the resulting system matrix \(\varPi A^{(\infty )}\varPi ^{\mathsf {T}}\) will be large scale and dense, which makes further computations barely feasible in terms of computation time and memory consumption. Also, a systematic error can easily be introduced by inaccurate computations with (20). Therefore, like in many similar applications, e.g., [2, 25, 28], we derive a suitable implicit implementation of the projection as it was proposed initially in [32].

We note that only the projected parts \(X_\varPi := \varPi ^{\mathsf {T}}X\varPi\) and \(Y_\varPi = \varPi ^{\mathsf {T}}Y\varPi\) of the Riccati solutions X and Y contribute to the controller; see (8) and [9]; where \(X_\varPi\) and \(Y_\varPi\) solve the equations

$$\begin{aligned}&\varPi \big (A^{(\infty )}\big )^{\mathsf {T}}\varPi ^{\mathsf {T}} X_\varPi E^{\mathsf {T}} + EX_\varPi \varPi A^{(\infty )}\varPi ^{\mathsf {T}} \nonumber \\ {}&\quad -(1 - \gamma ^{-2}) EX_\varPi B B^{\mathsf {T}} X_\varPi E^{\mathsf {T}} + \varPi C^{\mathsf {T}} C \varPi ^{\mathsf {T}} = 0, \end{aligned}$$
(23a)
$$\begin{aligned}&\varPi A^{(\infty )}\varPi ^{\mathsf {T}} Y_\varPi E + E^{\mathsf {T}}Y_\varPi \varPi \big (A^{(\infty )}\big )^{\mathsf {T}} \varPi ^{\mathsf {T}}\nonumber \\ {}&\quad - (1 - \gamma ^{-2}) E^{\mathsf {T}} Y_\varPi C^{\mathsf {T}} C Y_\varPi E + \varPi B B^{\mathsf {T}} \varPi ^{\mathsf {T}} = 0. \end{aligned}$$
(23b)

These Eqs. (23) are derived from (22) by pre- and postmultiplication with \(\varPi\) and \(\varPi ^{\mathsf {T}}\), respectively, and by means of the identities \(\varPi ^2 = \varPi\) and \(\varPi E=E\varPi ^{\mathsf {T}}\). In computational methods for large-scale sparse Riccati equations [12, 13, 34, 43, 44, 46], low-rank factors \(Z_k\) of the solution are computed such that, e.g., \(Y_\varPi \approx Z_k Z_k^{\mathsf {T}}\), by applying repeated solves of shifted linear systems that for (23b) read like

$$\begin{aligned} (\varPi A \varPi ^{\mathsf {T}} + p_i E) Z = W. \end{aligned}$$
(24)

With the requirement that the right-hand side lies in the appropriate subspace, i.e., \(W = \varPi ^{\mathsf {T}} W\), (24) can be equivalently formulated as

$$\begin{aligned} \left[ \begin{array}{ll} A + p_i E &{} J^{\mathsf {T}} \\ J &{} 0 \end{array}\right] \left[ \begin{array}{ll} Z \\ Z_\perp \end{array}\right]&= \left[ \begin{array}{ll} E^{\mathsf {T}} W \\ 0 \end{array}\right] . \end{aligned}$$

Here, the scalar \(p_i \in \mathbb {C}\) is a shift parameter that usually occurs in the considered iterative low-rank solvers like Newton-ADI or Krylov subspace methods; see, e.g., [12, 13, 34, 43, 44, 46]; and \(Z_\perp\) is an auxiliary variable.

3.2 Computation of low-rank controllers

The practical construction of a reduced-order controller for (18) follows, in principle, the different steps mentioned in this paper so far with some additional numerical tricks. For simplicity, we give a final summary of the performed steps in the following:

  • Step 1. Computation of a suitable robustness margin: So far, the margin \(\gamma\) was assumed to be given, but in fact it can be computed utilizing the necessary and sufficient conditions for the existence of a stabilizing controller from Proposition 1. In practice, we solve repeatedly (23) for different instances of the robustness margin and use (7) to determine the next iterate. Thereby, a sufficient \(\gamma\) can be computed.

  • Step 2. Construction of the reduced-order model: Now, the full-order system (18) needs to be reduced by Algorithm 1. Therefore, we use the final low-rank solution factors \(Z_{k}^{\mathsf {r}}\) and \(Z_{k}^{\mathsf {f}}\) of the projected algebraic Riccati equations (23a) and (23b), respectively, from the previous computation of the robustness margin in Step 1 corresponding to the computed \(\gamma\) such that \(X_\varPi \approx Z_{k}^{\mathsf {r}} \left( Z_{k}^{\mathsf {r}} \right) ^{\mathsf {T}}\) and \(Y_\varPi \approx Z_{k}^{\mathsf {f}} \left( Z_{k}^{\mathsf {f}} \right) ^{\mathsf {T}}\). Algorithm 1 is then used for the system matrices E, \(A^{(\infty )}\), B, C by setting the low-rank factors of the Riccati equations to be \(R = Z_{k}^{\mathsf {f}}\) and \(L = Z_{k}^{\mathsf {r}}\). The rest follows exactly Algorithm 1. The order r of the reduced-order model can be determined by using the criteria in Theorem 1 for the stabilization of the full-order system and possible disturbances. This gives the reduced-order system matrices \(\widehat{E}= I_{r}\), \(\widehat{A}\), \(\widehat{B}\) and \(\widehat{C}\).

  • Step 3. Construction of the reduced-order controller: The system matrices of the reduced-order central controller \(\widehat{K}\) can then be computed adapting the formulas in (8). First, an approximation to the solutions of the \(\mathcal {H}_{\infty }\)-Riccati equations (6) for the reduced-order system is directly given by

    $$\begin{aligned} \widehat{Y}_{\mathcal {H}_{\infty }} = W^{\mathsf {T}} E Z_{k}^{\mathsf {f}} (Z_{k}^{\mathsf {f}})^{\mathsf {T}} E^{\mathsf {T}} W \quad \text {and} \quad \widehat{X}_{\mathcal {H}_{\infty }} = T^{\mathsf {T}} E^{\mathsf {T}} Z_{k}^{\mathsf {r}} (Z_{k}^{\mathsf {r}})^{\mathsf {T}} E T, \end{aligned}$$

    with the low-rank factors \(Z_{k}^{\mathsf {f}}\) and \(Z_{k}^{\mathsf {r}}\) from the computation of the robustness margin, and W and T the truncation matrices from HINFBT in Algorithm 1. Then, the system matrices of the reduced-order controller are given by

    $$\begin{aligned} {\widetilde{\widehat{E}}}&= I_{r}, \end{aligned}$$
    (25a)
    $$\begin{aligned} {\widetilde{\widehat{A}}}&= \widehat{A}- (1 - \gamma ^{-2}) \widehat{Y}_{\mathcal {H}_{\infty }} \widehat{C}^{\mathsf {T}} \widehat{C}- \widehat{B}\widehat{B}^{\mathsf {T}} \widehat{X}_{\mathcal {H}_{\infty }} \widehat{Z}_{\mathcal {H}_{\infty }}, \end{aligned}$$
    (25b)
    $$\begin{aligned} {\widetilde{\widehat{B}}}&= \widehat{Y}_{\mathcal {H}_{\infty }} \widehat{C}^{\mathsf {T}}, \end{aligned}$$
    (25c)
    $$\begin{aligned} {\widetilde{\widehat{C}}}&= - \widehat{B}^{\mathsf {T}} \widehat{X}_{\mathcal {H}_{\infty }} \widehat{Z}_{\mathcal {H}_{\infty }}, \end{aligned}$$
    (25d)

    where \(\widehat{Z}_{\mathcal {H}_{\infty }} = (I_{r} - \gamma ^2 \widehat{Y}_{\mathcal {H}_{\infty }} \widehat{X}_{\mathcal {H}_{\infty }})^{-1}\).

4 Numerical examples

4.1 Example setups

As numerical examples, we consider here two-dimensional flows through a channel with circular obstacles, exemplarily depicted in Fig. 1, with controls acting on the boundary of the obstacles, and with observation of locally spatially averaged velocities in a domain of observation downstream of the obstacles. The generic model then reads

$$\begin{aligned} \dot{v}(t) +(v(t)\cdot \nabla )v(t) - \frac{1}{\mathsf {Re}} \varDelta v(t) + \nabla p(t)&= 0, \quad \text {in}\quad \varOmega , \end{aligned}$$
(26a)
$$\begin{aligned} \nabla \cdot v(t)&= 0, \quad \text {in}\quad \varOmega , \end{aligned}$$
(26b)
$$\begin{aligned} \tfrac{1}{\mathsf {Re}}\tfrac{\partial v}{\partial n}v(t) -np(t)&= 0, \quad \text {on}\quad \varGamma _{\text {out}}, \end{aligned}$$
(26c)
$$\begin{aligned} v(t)&= 0, \quad \text {on}\quad \varGamma _w, \end{aligned}$$
(26d)
$$\begin{aligned} v(t)&= ng_{\text {in}}, \quad \text {on}\quad \varGamma _{\text {in}}, \end{aligned}$$
(26e)
$$\begin{aligned} v(t)&= ng_iu_i(t), \quad \text {on}\quad \varGamma _i, \quad i=1,2. \end{aligned}$$
(26f)

Here, n denotes the inward normal of the boundaries, \(g_{\text {in}}\) models the inflow boundary condition via a parabola that takes the value 0 at the edges and 2/3 at the center of the inflow boundary, respectively, and \(g_i\) and \(u_i\) are control shape functions and scalar control value functions. The condition (26c) for the outflow is the standard do-nothing condition. The geometrical extensions, the choice of the parameters \(u_{\text {in}}\) and \(\mathsf {Re}\), and the definition of the output operator \(\mathcal {C}\), and of the shape functions are given in the description of the test cases below and in Table 1.

Fig. 1
figure 1

Computational domain of the cylinder wake

Full size image

4.1.1 Test case: cylinder wake

As first example, we consider the cylinder as it has been described in the benchmark work [42] in an enlarged domain to reduce the stabilizing effects of the channel walls; see Fig. 2. To readily include the boundary controls in the finite elements discretization, we relax the Dirichlet control conditions (26f) towards Neumann conditions via

$$\begin{aligned} v(t) = ng_i u_i(t) - \alpha \left( \tfrac{1}{\mathsf {Re}}\tfrac{\partial v}{\partial n}(t) -np(t)\right) \quad \text {on}\quad \Gamma _i, \quad i=1,2, \end{aligned}$$
(27)

with a parameter \(\alpha\) that is supposed to be small; see, e.g., [33] for convergence properties of this relaxation in optimal control of stationary flows. The shape functions \(g_i\) are chosen to be \(g_i(x) = \cos (s_i(x)) - 1\), where \(s_i:\Gamma _i \rightarrow [0, 2\pi ]\), for \(i=1,2\), parametrizes the arc length of the control boundaries; see [5, Sec. 9.3], where also the assembling of the associated control operator is explained.

Fig. 2
figure 2

Example simulation of the cylinder wake with no control input at \(\mathsf {Re}= 60\)

Full size image

4.1.2 Test case: double cylinder

As second example, we borrow the numerical setup from [18] of the wake with two cylinders in two dimensions; see Fig. 3. The actuation of the flow happens through controlled rotation of the individual cylinders. This means that, instead of the inflow conditions (26f), we prescribe the control boundary conditions as

$$\begin{aligned} v\bigr |_{\Gamma _i} = t_{\Gamma _i}r_i u_i , \quad i=1,2, \end{aligned}$$

where the \(\Gamma _i\) are the two boundaries of the cylinders, \(t_{\varGamma _i}\) denote the tangential vectors at the boundaries, \(r_i\) are the radii, and the \(u_i\) are the two scalar input functions, for \(i=1,2\). Although these boundary conditions have no normal component, e.g., \(v\cdot n = 0\) at \(\varGamma _i\), for \(i=1,2\), and, thus, can be included in a weak formulation with a bounded input operator \(\mathcal {B}\) (cf. [4, Sec. 3.1]), we use the same Robin relaxation (5, Sec. 9.3] with the penalization parameter 27) that is defined to relax controls with normal components.

Fig. 3
figure 3

Example simulation of the double cylinder example with no control input at \(\mathsf {Re}= 60\)

Full size image

4.1.3 Further computational setup for both test cases

The spatial discretization in both test cases is done with \(P_2-P_1\) Taylor-Hood finite elements. The input operators are assembled as described in [\(\alpha\) set to \(10^{-5}\). Since the output operator \(\mathcal {C}\) is of distributed type, its assembling with finite elements is straight forward.

We assemble the coefficients of the corresponding linearized system (19) and follow the procedure laid out in Sect. 3.2. Therefor, we compute robustness margins \(\gamma\) and the corresponding low-rank approximations to the stabilizing solutions of (23) using the low-rank Riccati iteration method from [15]. This approach is suited to solve Riccati equations like (6) with indefinite quadratic terms by splitting the computations into two steps: first, the solution to the classical LQG-Riccati equations is computed and, afterwards, residual Riccati equations are solved to update the overall solution. While this method was originally developed for Riccati equations with indefinite quadratic terms, it also allows us here to efficiently compute the stabilizing solutions of (23) for many different instances of the robustness margin to find suitable \(\gamma\) for controller design. Therein, we use the low-rank Newton-ADI [12, 46] as solver for the LQG-Riccati equations and the low-rank RADI [13] for the update residual equations. The final computed robustness margins \(\gamma\) for both test examples are given in Table 2.

The reduced central output-based feedback controller is then synthesized via (25), with its coefficients

$$\begin{aligned} \widetilde{\widehat{A}}, \; \widetilde{\widehat{B}}, \;\text {and}\;\widetilde{\widehat{C}}. \end{aligned}$$

Thus, with \(v_\infty\) denoting the target state and linearization point, the closed-loop system reads

$$\begin{aligned} E \dot{v}&= A_{\mathsf {S}}v + F(v) + J^{\mathsf {T}}p + B\widetilde{\widehat{C}}\widetilde{\widehat{x}}+f, \\ 0&= Jv + g, \\ \dot{\widetilde{\widehat{x}}}&= \widetilde{\widehat{A}}\widetilde{\widehat{x}}+ \widetilde{\widehat{B}}C(v - v_\infty ). \end{aligned}$$

For the time discretization, we use a uniform grid of size h and employ backward differencing of second order in the linear part including the controller and the extrapolation

$$\begin{aligned} F(v(t+h)) \approx 2F(v(t)) - F(v(t-h)) \end{aligned}$$

for the nonlinearity. Using one step of Heun’s method for the initialization, this results in a time integration scheme of order two. As the initial value, we use the corresponding steady state \(v_\infty\). To trigger the instabilities, we add an input perturbation that acts at the beginning of the simulation like

$$\begin{aligned} u_{\delta } = {\left\{ \begin{array}{ll} 10^{-6}\sin (2\pi t), &{} \quad \text {for}\;t \in [0,1], \\ 0, &{} \quad \text {for}\;t > 1. \end{array}\right. } \end{aligned}$$

The final simulation setups and all parameters that define the simulations for both test cases are listed in Tables 1 and 2.

Table 1 Simulation setups
Full size table
Table 2 Simulation and controller parameters
Full size table

Remark 3

(Validity of Assumption 1) The problem at hand derives from a state-space system with input operator B and output operator C that is brought into the normalized form (9). Accordingly, Assumption 1 is reduced to stability and detectability with respect to the given inputs and outputs; cf. the discussion after Eq. (9). Due to the involvement of the projector \(\varPi\) in the underlying ODE (21), for example, the system \((sE-\varPi A^{(\infty )}\varPi ^{\mathsf {T}}, \varPi B)\) will always have nonstabilizable modes on the imaginary axis so that it can not be stabilizable in the standard definition. However, the modes associated with the kernel of \(\varPi\) are excluded from the dynamics, so that notions of stabilizability and detectability can be adapted, e.g., by using the reduced coordinates as in [32].

Still, an analytical confirmation of Assumption 1 is a difficult task. There exist relevant fundamental work (see, e.g. [38]) but these results are generic and do not respect particular input operators and specific domains, let alone the discretization. Instead, one may resort to a numerical approach to establish stabilizability as proposed in [7, Sec. 5].

Finally, we want to remark that this stability and detectability is both necessary and sufficient for the existence of the stabilizing Riccati solutions and the convergence of the employed algorithms. Thus, a failure of the algorithms means that the assumptions do not hold and vice versa.

4.2 Numerical results

The experiments reported here have been executed on a machine with 2 Intel(R) Xeon(R) Silver 4110 CPU processors running at 2.10 GHz and equipped with 192 GB total main memory. The computer is run on CentOS Linux release 7.5.1804 (Core).

The solutions to the large-scale projected Riccati equations (23) have been computed in MATLAB 9.7.0.1190202 (R2019b) using the routines from the M-M.E.S.S. library version 2.0.1 [16, 41]. Also, the errors in coprime factorizations have been computed in MATLAB. For the spatial discretization, we employ the FEniCS [1] finite elements toolbox and the python module dolfin_navier_scipy [29] to extract the discrete operators for the computation of the feedback gains and for the time integration, which is done in SciPy.

We investigate the computed controllers in terms of robustness against controller reduction and robustness against linearization errors via the following criteria:

  • Does the robustness margin \(\gamma\) cover the linearization error, i.e.,

    $$\begin{aligned} \left\Vert \varDelta \right\Vert _{\mathcal {H}_{\infty }} := \left\Vert \left[ \begin{array}{ll} N - N_{\varDelta }&\; M - M_{\varDelta } \end{array}\right] \right\Vert _{\mathcal {H}_{\infty }}&< \gamma ^{-1} \end{aligned}$$

    as in (17)?

  • Does the robustness margin cover the a-priori estimate for the truncation error, i.e., \(\epsilon ( \beta + \gamma ) < 1\) as in Theorem 1? (For reference, we also check the less conservative a-posteriori estimate \(\hat{\epsilon }(\beta + \gamma ) < 1\); cf. (13).)

  • Does the controller work in the numerical experiment, i.e., does it stabilize the steady state by completely suppressing the oscillations that can be observed in the uncontrolled cases as illustrated in Figs. 2 and 3?

Remark 4

(Comments on computing \(\epsilon\)) In the presented setup where the solutions to the Riccati equations are approximated via low-rank factorizations, the value of most of the discarded characteristic \(\mathcal {H}_{\infty }\) values \(\sigma _k\) is not known so that the precise computation of \(\epsilon\) as in (14) is not possible. As it is common practise, we set those \(\sigma _k\) that are not covered from the low-rank factorization to zero and compute \(\epsilon\) by means of the remaining. This is well inline with the working hypothesis of the low-rank approximations, that the characteristic values show a rapid decay towards zero.

Remark 5

(Comments on linearization errors) Our analysis of the linearization error follows the goal of robustifying the observer-based controller. Thus, we will consider linearization errors like inexact computations of the linearization points as it may arise in numerical calculations. Whether or not the linear controller is suitable for the nonlinear model is a question of performance, which is not considered here. In our numerical examples, we choose the linearization point as starting point and add small input perturbations such that we may well assume that we are close to the working point and the linearization principle for controller design ensures the performance. For estimates on the admissible deviation from the working point in the infinite-dimensional Navier–Stokes model see, e.g., [40].

To examine the linearization error, we proceed as follows: For the cylinder wake, we consider the linearization error that stems from an incomplete iteration for the steady state computation. That is, in (21), instead of the exact linearization \(A^{(\infty )}=A_{\mathsf {S}}+ (\partial _v F)(v_\infty )\) based on the exact steady state \(v_\infty\), we consider \(A^{(\ell )}:= A_{\mathsf {S}}+ (\partial _v F)(v_\ell )\), where \(v_\ell\) is the approximation to \(v_\infty\) that is obtained after \(\ell\) Picard iterations started from the steady state solution for \(\mathsf {Re}=40\). For the double cylinder, we found that at \(\mathsf {Re}= 60\), the steady state is so unstable that the Picard iteration does not converge. The computation can be done by a Newton iteration that, however, does not provide a smooth parametrization of the approximation in the relevant region because of its fast convergence. Therefor, we consider the perturbed coefficient \(A^{(\ell )}= A_{\mathsf {S}}+ (\partial _v F)(v^{(\ell )})\), where \(v^{(\ell )}\) as the steady state solution to the problem with the Reynolds number perturbed by \(\ell\) thousandths, i.e., \(\mathsf {Re}^{(\ell )} = (1 + \frac{\ell }{1000}) \cdot 60\). We parametrize the truncation error via the truncation threshold tol that defines \(\texttt {tol}> \sigma _{r+1} \ge \sigma _{r+2} \ge ...\), i.e., the size of those characteristic \(\mathcal {H}_{\infty }\)-values that are discarded in Algorithm 1.

We have checked the performance in the simulations by examining the empirical variances in the time series of the output signal in the third and fourth quarters of the time interval. If the difference between the fourth and third segment is negative, we conclude that the oscillations were on the decline. If the difference is positive but in the order of \(10^{-15}\), we conclude that the signals were dominated by numerical errors. In both cases, the corresponding setup was reported as stabilizing in the simulation. In Table 3, we have tabulated the computed coprime factor errors of the different linearizations and identified the threshold where \(\Vert \varDelta _\ell \Vert < \gamma ^{-1}\).

Table 3 Computed left coprime factor errors
Full size table

The search of the parameter ranges for successful controller setups is illustrated in Fig. 4. Therein, each mark denotes a simulation that either failed (unstable behavior) or has been stabilized. For a better practical understanding, we added additionally to the truncation tolerance tol the orders of the reduced controllers used in the simulations. The a-priori and a-posteriori stability criteria are added as horizontal lines and the bound for the linearization error as vertical line. The intersection of both hatched areas are the regions where the controller setups are predicted to be stable by both criteria.

Fig. 4
figure 4

Performance of reduced controllers for the test examples

Full size image

For the cylinder wake, the estimates precisely confine the range where the controllers are functional; see Fig. 4a. The a-priori estimate (15) turns out to be away from the observed failures by a factor larger than 8, and the bound on the linearization error by a factor of \(2^{\frac{1}{2}}\). Successful stabilization has been observed further outside of the region predicted by the bound on the linearization error and also outside the a-posteriori estimate defined by the truncation error.

The simulations for the double cylinder are shown in Fig. 4b. Here, the picture is even more compliant with the predictions as the a-posteriori estimate for the region of stabilization perfectly separates the failed from the stabilized simulation. The a-priori bound for the truncation error suggests a tol of 0.002, which is a significant underestimate of the stability region and therefor not shown in Fig. 4b. On the other hand, the bound for the linearization is again too conservative for the actual region of stabilization, this time by a factor of 4. Curiously, increasing slightly the order of the \(\mathcal {H}_{\infty }\)-controller allows to perform successful simulations with a perturbation of the \(\mathsf {Re}\) number up to about 14%. As for the cylinder wake, this shows that successful stabilization can be observed outside the region predicted by the linearization error bound while the error estimates for the truncation error do not leave a margin.

5 Conclusions

As the numerical examples have illustrated, the use of the general Riccati-based low-order \(\mathcal {H}_{\infty }\)-controller design is well feasible for incompressible flows in simulations. The provided estimates on the guaranteed robustness have been proven reliable though, in some cases, conservative. Together with the theoretical results of earlier works that \(\mathcal {H}_{\infty }\)-controller can compensate various model errors, the availability of efficient general purpose numerical methods is key for the applicability of these model-based controllers in simulations and even experiments.