2.1. System definition

We consider the linear time-invariant system

(2.1)\begin{equation} \left.\begin{array}{c@{}}\dfrac{{\rm d} \boldsymbol{u}}{{\rm d} t} (t) ={{\boldsymbol{\mathsf{A}}}} \boldsymbol{u}(t) +{{{\boldsymbol{\mathsf{B}}}}_{\boldsymbol{f}}}\boldsymbol{f}(t) +{{{\boldsymbol{\mathsf{B}}}}_{\boldsymbol{a}}}\boldsymbol{a}(t), \\ \boldsymbol{y}(t) ={{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol{y}}}\boldsymbol{u}(t) + \boldsymbol{n}(t),\\ \boldsymbol{z}(t)= {{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol{z}}}\boldsymbol{u}(t), \end{array}\right\}\end{equation}

where $\boldsymbol {u} \in \mathbb {C}^{n_u}$ represents the flow state, $\boldsymbol {f} \in \mathbb {C}^{n_f}$ is an unknown stochastic forcing, which can represent external disturbances and/or nonlinear interactions (McKeon & Sharma Reference McKeon and Sharma2010), $\boldsymbol {a}\in \mathbb {C}^{n_a}$ represents flow actuation used for control, $\boldsymbol {y} \in \mathbb {C}^{n_y}$ is a set of system observables, $\boldsymbol {n} \in \mathbb {C}^{n_y}$ is measurement noise and $\boldsymbol {z} \in \mathbb {C}^{n_z}$ is a set of targets for the control problem, for instance, perturbations at a given position or surface loads, to be minimized. The system evolution is described by the matrix ${{\boldsymbol{\mathsf{A}}}} \in \mathbb {C}^{n_u\times n_u}$, representing the linearized Navier–Stokes operator. The matrices ${{{\boldsymbol{\mathsf{B}}}}_{\boldsymbol {f}}} \in \mathbb {C}^{n_u\times n_f}, {{{\boldsymbol{\mathsf{B}}}}_{\boldsymbol {a}}} \in \mathbb {C}^{n_u\times n_a}, {{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol {y}}} \in \mathbb {C}^{n_y\times n_u}$ and ${{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol {z}}} \in \mathbb {C}^{n_z\times n_u}$ determine the spatial support of external disturbances, actuators, sensors and targets, as in Bagheri, Brandt & Henningson (Reference Bagheri, Brandt and Henningson2009a). Forcing and sensor noise are modelled as stochastic zero-mean processes with two-point space–time correlations given by

(2.2)$$\begin{gather} {{\boldsymbol{\mathsf{F}}}}(t-t') = \langle \boldsymbol{f}(t)\boldsymbol{f}^{{{\dagger}}}(t') \rangle , \end{gather}$$

(2.3)$$\begin{gather}{{\boldsymbol{\mathsf{N}}}}(t-t') = \langle \boldsymbol{n}(t)\boldsymbol{n}^{{{\dagger}}} (t')\rangle, \end{gather}$$

with ${\dagger}$ representing the adjoint operator using a suitable inner product and $\langle \cdot \rangle$ representing the ensemble average. We emphasize that these forcing and noise statistics are more general than those assumed in the derivation of the Kalman filter and LQG control, in which they must be uncorrelated in time. A frequency-domain representation of the two-point correlation is given by the CSDs

(2.4)$$\begin{gather} \hat{{{\boldsymbol{\mathsf{F}}}}}(\omega) = \langle \hat{\boldsymbol{f}} \hat{\boldsymbol{f}}^{{{\dagger}}} \rangle, \end{gather}$$

(2.5)$$\begin{gather}\hat{\boldsymbol{\mathsf{N}}} (\omega) = \langle \hat{\boldsymbol{n}} \hat{\boldsymbol{n}}^{\dagger} \rangle. \end{gather}$$

Note that this differs from the standard definition, i.e. $\langle \hat {\boldsymbol {f}} \hat {\boldsymbol {f}}^{H} \rangle$, with $H$ representing the conjugate transpose of the matrices; this modified definition simplifies the derivations that follow.

Throughout this work, we assume the system to be stable, that is, all eigenvalues $\lambda$ of ${{\boldsymbol{\mathsf{A}}}}$, i.e.  $\lambda$ for which ${{\boldsymbol{\mathsf{A}}}} \hat {\boldsymbol {u}} = -\mathrm {i} \lambda \hat {\boldsymbol {u}}$ has non-trivial solutions, lie in the lower half-plane. As discussed in § 1, optimal control strategies are best suited for amplifier flows, which satisfy this stability requirement.

It is useful to split (2.1) into two systems, one of which is driven by the forcing and includes sensor noise,

(2.6) \begin{equation} \left.\begin{array}{c@{}}\dfrac{{\rm d} \boldsymbol{u}_1}{{\rm d} t}(t) ={{\boldsymbol{\mathsf{A}}}} \boldsymbol{u}_1(t) +{{{\boldsymbol{\mathsf{B}}}}_{\boldsymbol{f}}}\boldsymbol{f}(t), \\ \boldsymbol{y}_1(t) ={{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol{y}}}\boldsymbol{u}_1(t) + \boldsymbol{n}(t),\\ \boldsymbol{z}_1(t) = {{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol{z}}} \boldsymbol{u}_1(t), \end{array}\right\}\end{equation}

and another noiseless system driven by actuation only,

(2.7)\begin{equation} \left.\begin{array}{c@{}} \dfrac{{\rm d} \boldsymbol{u}_2}{{\rm d} t}(t) = {{\boldsymbol{\mathsf{A}}}} \boldsymbol{u}_2(t) + {{{\boldsymbol{\mathsf{B}}}}_{\boldsymbol{a}}} \boldsymbol{a}(t), \\ \boldsymbol{y}_2(t) = {{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol{y}}} \boldsymbol{u}_2(t) ,\\ \boldsymbol{z}_2(t) = {{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol{z}}} \boldsymbol{u}_2(t). \end{array}\right\} \end{equation}

The original system can be recovered adding the variables with subscripts ‘1’ and ‘2’, that is,

(2.8a–c)\begin{equation} \boldsymbol{u} = \boldsymbol{u}_1 + \boldsymbol{u}_2, \quad \boldsymbol{y} = \boldsymbol{y}_1 + \boldsymbol{y}_2, \quad \boldsymbol{z} = \boldsymbol{z}_1 + \boldsymbol{z}_2. \end{equation}

Figures 1 and 2 illustrate these systems.

Figure 1.

Block diagram representing the system (2.1).

Figure 2.

Block diagram representing the system (2.6) and (2.7).

The derivations that follow make use of the frequency-domain form of these equations. Taking a Fourier transform of (2.6) and (2.7) yields the input–output relationships

(2.9)$$\begin{gather} \hat{\boldsymbol{y}}_1(\omega) = {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}\boldsymbol{f}}(\omega) \hat{\boldsymbol{f}}(\omega) + \hat{\boldsymbol{n}}(\omega) , \quad \text{with } {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}\boldsymbol{f}}(\omega) = {{{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol{y}}} {{\boldsymbol{\mathsf{R}}}} {{{\boldsymbol{\mathsf{B}}}}_{\boldsymbol{f}}}}, \end{gather}$$

(2.10)$$\begin{gather}\hat{\boldsymbol{z}}_1(\omega) = {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{z}\boldsymbol{f}}(\omega)\hat{\boldsymbol{f}}(\omega), \quad \text{with } {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{z}\boldsymbol{f}}(\omega) = {{{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol{z}}} {{\boldsymbol{\mathsf{R}}}} {{{\boldsymbol{\mathsf{B}}}}_{\boldsymbol{f}}}}, \end{gather}$$

(2.11)$$\begin{gather}\hat{\boldsymbol{y}}_2(\omega) = {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}\boldsymbol{a}}(\omega) \hat{\boldsymbol{a}}(\omega), \quad \text{with } {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}\boldsymbol{a}}(\omega) = {{{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol{y}}} {{\boldsymbol{\mathsf{R}}}} {{{\boldsymbol{\mathsf{B}}}}_{\boldsymbol{a}}}}, \end{gather}$$

(2.12)$$\begin{gather}\hat{\boldsymbol{z}}_2(\omega) = {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{z}\boldsymbol{a}}(\omega)\hat{\boldsymbol{a}}(\omega), \quad \text{with } {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{z}\boldsymbol{a}}(\omega) = {{{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol{z}}} {{\boldsymbol{\mathsf{R}}}} {{{\boldsymbol{\mathsf{B}}}}_{\boldsymbol{a}}}}, \end{gather}$$

where

(2.13)\begin{equation} {{\boldsymbol{\mathsf{R}}}} = (-\mathrm{i}\omega{{\boldsymbol{\mathsf{I}}}}-{{\boldsymbol{\mathsf{A}}}})^{{-}1} \end{equation}

is the resolvent operator. Note that the above input–output relations are exact within the linearized model used here, and thus has the same limitation as any linear modelling of the original nonlinear system, i.e. treating the nonlinear interactions as exogenous stochastic noise. This framework is shared by all linear control laws that have been developed.

2.2. The estimation problem

Before deriving the optimal causal estimator that we seek, we first briefly review the derivation of the optimal non-causal estimator (Martini et al. Reference Martini, Cavalieri, Jordan, Towne and Lesshafft2020). We will see later that the two cases are closely related but that imposing causality leads to the appearance of an additional term. We focus on the uncontrolled estimation problem, i.e. $\boldsymbol {a}(t) = 0$, thus making (2.1) and (2.6) equivalent.

The optimal estimator minimizes the cost functional

(2.14)\begin{equation} J = \int_{-\infty}^{\infty} \langle \boldsymbol{e}^{{{\dagger}}}(t) \boldsymbol{e}(t) \rangle \, {\rm d} t = \frac{1}{2{\rm \pi}} \int_{-\infty}^{\infty} \langle \hat{\boldsymbol{e}}^{{{\dagger}}}(\omega) \hat{\boldsymbol{e}}(\omega) \rangle \, {\rm d} \omega, \end{equation}

defined in terms of the estimation error

(2.15)\begin{equation} \boldsymbol{e}(t) = \boldsymbol{z}(t) - \tilde{\boldsymbol{z}}(t), \end{equation}

where $\boldsymbol {z}$ is the target, i.e. the quantity to be estimated, and $\tilde {\boldsymbol {z}}$ is its estimate and where $\langle\, \cdot\, \rangle$ represents the expected value. Note that full-state estimation is recovered using ${{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol {z}}}={{\boldsymbol{\mathsf{I}}}}$.

We seek an estimate of the targets in the form of a linear combination of sensor readings, in the form of

(2.16a,b)\begin{equation} \tilde{\boldsymbol{z}}(t) = \int_{-\infty}^{\infty}\boldsymbol{\mathsf{T}}_{z,nc}(\tau) \boldsymbol{y}(t-\tau) \, {\rm d} \tau,\quad \tilde{\hat{\boldsymbol{z}}} (\omega) = \hat{\boldsymbol{\mathsf{T}}}_{z,nc} \boldsymbol{\hat{y}}(\omega), \end{equation}

where the subscript $nc$ is a reminder that the kernel, ${{\boldsymbol{\mathsf{T}}}}_{z,nc} \in \mathbb {C}^{n_z\times n_y }$, is, in general, non-causal, and $\hat {{{\boldsymbol{\mathsf{T}}}}}_{z,nc}$ is its Fourier transform. Expanding (2.14) leads to

(2.17)\begin{equation} J = \frac{1}{2{\rm \pi}} \int_{-\infty}^{\infty} (( {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{z}\boldsymbol{f}} - \hat{{{\boldsymbol{\mathsf{T}}}}}_{z,nc} {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}\boldsymbol{f}} )^{{{\dagger}}} \hat{{{\boldsymbol{\mathsf{F}}}}} ({{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{z}\boldsymbol{f}} - \hat{{{\boldsymbol{\mathsf{T}}}}}_{z,nc} {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}\boldsymbol{f}} ) + {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}\boldsymbol{f}}^{{{\dagger}}} \hat{{{\boldsymbol{\mathsf{T}}}}}_{z,nc}^{{{\dagger}}} \hat{{{\boldsymbol{\mathsf{N}}}}}\hat{{{\boldsymbol{\mathsf{T}}}}}_{z,nc} {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}\boldsymbol{f}} ) \, {\rm d} \omega. \end{equation}

In the above equations, and henceforth, the frequency dependence is omitted for clarity. The minimum is found by taking the derivative of the cost function (2.17) with respect to $\hat {{{\boldsymbol{\mathsf{T}}}}}_{z,nc}^{{\dagger} }(\omega )$ and setting it to zero. The optimal estimation kernel is obtained as the solution of

(2.18)\begin{equation} \hat{{{\boldsymbol{\mathsf{T}}}}}_{z,nc} \hat{{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{l}}}}}} = \hat{{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{r}}}}}}, \end{equation}

where

(2.19)$$\begin{gather} \hat{{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{l}}}}}} = {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}\boldsymbol{f}} \hat{{{\boldsymbol{\mathsf{F}}}}} {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}\boldsymbol{f}}^{{{\dagger}}} + \hat{\boldsymbol{\mathsf{N}}}, \end{gather}$$

(2.20)$$\begin{gather}\hat{{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{r}}}}}}= {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{z}\boldsymbol{f}} \hat{{{\boldsymbol{\mathsf{F}}}}} {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}\boldsymbol{f}} ^{{{\dagger}}}. \end{gather}$$

Note that the forcing CSD $\hat {{{\boldsymbol{\mathsf{F}}}}}$ appears explicitly in the equation, and is thus naturally handled by the approach.

Causality of the estimation kernel can be enforced in (2.14) using Lagrange multipliers (Martinelli Reference Martinelli2009), as

(2.21)\begin{align} J' &= J + \int_{-\infty}^{\infty} {\rm Tr}( \boldsymbol{\varLambda}_-(t) {{\boldsymbol{\mathsf{T}}}}_{z,c}(t)+\boldsymbol{\varLambda}^{{{\dagger}}}_- (t) {{\boldsymbol{\mathsf{T}}}}^{{{\dagger}}}_{z,c}(t))\,{\rm d} t \nonumber\\ &= J + \int_{-\infty}^{\infty} {\rm Tr}(\hat{\boldsymbol{\varLambda}}_- \hat{{{\boldsymbol{\mathsf{T}}}}}_{z,c}+ \hat{\boldsymbol{\varLambda}}^{{{\dagger}}}_- \hat{{{\boldsymbol{\mathsf{T}}}}} ^{{{\dagger}}}_{z,c}) \,{\rm d} \omega, \end{align}

where the Lagrange multipliers ${\boldsymbol {\varLambda }}_-,{\boldsymbol {\varLambda }}_+ \in \mathbb {C}^{n_y\times n_z}$ are required to be zero for $t > 0$, thus enforcing the condition ${{{\boldsymbol{\mathsf{T}}}}_{z,c}(t<0)=0 }$. The subscript $c$ is used to emphasize the causal nature of this kernel.

The functionals $J$ and $J'$ differ only by linear terms in the estimation kernel, and thus it is straightforward to see that taking the derivative of $J'$ and setting it to zero leads to

(2.22)\begin{equation} \hat{{{\boldsymbol{\mathsf{T}}}}}_{z,c} \hat{{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{l}}}}}} + \hat{\boldsymbol{\varLambda}}_- = \hat{{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{r}}}}}}. \end{equation}

This apparently simple equation hides significant complexities. First, this is a single equation with two variables, ${\hat {{{\boldsymbol{\mathsf{T}}}}}_{z,c}}$ and ${\hat {\boldsymbol {\varLambda }}}$. Nevertheless, it admits a unique solution due to the requirements that ${{{\boldsymbol{\mathsf{T}}}}_{z,c}(t<0)=0 }$ and ${\boldsymbol {\varLambda }(t>0)=0 }$, which in the frequency domain impose restrictions on ${\hat {{{\boldsymbol{\mathsf{T}}}}}_{z,c}}$ and ${\hat {\boldsymbol {\varLambda }}}$, namely that these quantities are regular on the upper and lower complex planes, respectively. This restriction means that the values of ${\hat {{{\boldsymbol{\mathsf{T}}}}}_{z,c}}$ and ${\hat {\boldsymbol {\varLambda }}}$ for different frequencies have a non-trivial relation between them, and thus cannot be chosen independently. Equation (2.22), with the regularity constrains, constitutes a Wiener–Hopf problem. As discussed in § 1, analytical solutions are only known for special cases, and thus we resort to numerical methods. An introduction to Wiener–Hopf problems and numerical methods to solve them is presented in Appendix A. The solution of (2.22), once inverse Fourier transformed, is the optimal causal estimation kernel for the linear system at hand.

In previous studies (Martini et al. Reference Martini, Cavalieri, Jordan, Towne and Lesshafft2020; Amaral et al. Reference Amaral, Cavalieri, Martini, Jordan and Towne2021) we have shown that using the spatiotemporal forcing statistics considerably improves the accuracy of the estimation of a turbulent channel flow. The estimation method presented here preserves the ability to handle these complex forcing models, while being applicable in real time via the simple integration of

(2.23)\begin{equation} \tilde{\boldsymbol{z}}(t) = \int_{0}^{\infty} \boldsymbol{\mathsf{T}}_{z,c}(\tau) \boldsymbol{y}(t-\tau) \, {\rm d} \tau, \end{equation}

which is similar to (2.16a,b), but with integration restricted to positive $\tau$, implying that only present and past sensor measurements are used for estimation.

2.3. Partial-knowledge control

Control can be divided in full-knowledge control, where the full state of the flow is assumed to be known, and partial-knowledge control, where the flow state needs to be estimated from limited, noisy, sensor readings. In this section we derive the optimal partial-knowledge control; the full-knowledge case is considered in Appendix B for completeness.

Analogous to the estimation problem, we seek a control law that constructs an actuation signal as a linear function of sensor readings,

(2.24)\begin{equation} \boldsymbol{a}(t) = \int_{-\infty}^{\infty} {\boldsymbol{\varGamma}}'(\tau) {\boldsymbol{y}} (t-\tau) \, {\rm d} \tau, \quad \hat{\boldsymbol{a}} = \hat{\boldsymbol{\varGamma}}' \hat{\boldsymbol{y}}, \end{equation}

where $\boldsymbol {\varGamma }' \in \mathbb {C}^{ n_a\times n_y }$ is the control kernel and $\hat {\boldsymbol {\varGamma }}'$ is its Fourier transform. However, the derivation of the kernel is simplified if instead the actuation is expressed only in terms of $\boldsymbol {y}_1$, i.e. ignoring the influence of the actuator on the sensor,

(2.25)\begin{equation} \boldsymbol{a}(t) = \int_{-\infty}^{\infty} {\boldsymbol{\varGamma}}(\tau) {\boldsymbol{y}}_1 (t-\tau) \, {\rm d} \tau, \quad \hat{\boldsymbol{a}} = \hat{\boldsymbol{\boldsymbol{\varGamma}}} \hat{\boldsymbol{y}}_1 , \end{equation}

where again $\boldsymbol {\varGamma } \in \mathbb {C}^{ n_a\times n_y }$ and $\hat {\boldsymbol {\varGamma }}$ its Fourier transform.

This implies no loss of generality: as $\boldsymbol {y}_2$ is a function only of the previous actuation, which is known, it can be computed using the actuator impulse response and subtracted from $\boldsymbol {y}$ to obtain $\boldsymbol {y}_1$. The approach is equivalent to the formalism of an internal model control (Morari & Zafiriou Reference Morari and Zafiriou1989). The control kernel $\boldsymbol {\varGamma }$ is then chosen so as to minimize a cost functional that trades off the expected value of targets and actuation,

(2.26)\begin{equation} J = \int_{-\infty}^{\infty} \langle \boldsymbol{z}^{{{\dagger}}} (t) \boldsymbol{z}(t) + \boldsymbol{a}^{{{\dagger}}}(t){{\boldsymbol{\mathsf{P}}}}\boldsymbol{a}(t) \rangle \, {\rm d} t = \int_{-\infty}^{\infty} \langle \hat{\boldsymbol{z}}^{{{\dagger}}} \hat{\boldsymbol{z}} + \hat{\boldsymbol{a}}^{{{\dagger}}}{{\boldsymbol{\mathsf{P}}}}\hat{\boldsymbol{a}} \rangle \, {\rm d} \omega, \end{equation}

where ${{\boldsymbol{\mathsf{P}}}} \in \mathbb {C}^{ n_a\times n_a }$ is a positive-definite matrix containing actuation penalties. For simplicity, we do not include a state cost matrix. This does not imply any loss of generality, as any state cost matrix ${{\boldsymbol{\mathsf{Q}}}} \in \mathbb {C}^{ n_u\times n_u }$, which must be positive definite, can be absorbed into the definition of the targets as ${{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol {z}}}'= {{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol {z}}} {{\boldsymbol{\mathsf{Q}}}}_c$, where ${{\boldsymbol{\mathsf{Q}}}} = {{\boldsymbol{\mathsf{Q}}}}_c {{\boldsymbol{\mathsf{Q}}}}_c^{H}$ is a Cholesky decomposition of ${{\boldsymbol{\mathsf{Q}}}}$.

Using the identity $\textrm {Tr}(\boldsymbol {\varPhi }\boldsymbol {\varPsi })=\textrm {Tr}(\boldsymbol {\varPsi }\boldsymbol {\varPhi } )$, valid for any matrices $\boldsymbol {\varPsi }$ and $\boldsymbol {\varPhi }$ with suitable sizes, the functional is rewritten as

(2.27)\begin{equation} J =\int_{-\infty}^{\infty} \langle {\rm Tr}( \hat{\boldsymbol{z}}(\omega)\hat{\boldsymbol{z}}^{{{\dagger}}} (\omega)) + {\rm Tr}({{\boldsymbol{\mathsf{P}}}}\hat{\boldsymbol{a}}(\omega)\hat{\boldsymbol{a}}^{{{\dagger}}}(\omega)) \rangle \, {\rm d} \omega. \end{equation}

The functional is expanded using

(2.28)$$\begin{gather} \langle \hat{\boldsymbol{z}} \hat{\boldsymbol{z}}^{{{\dagger}}} \rangle = \langle \hat{\boldsymbol{z}}_1 \hat{\boldsymbol{z}}_1^{{{\dagger}}} \rangle+ \langle \hat{\boldsymbol{z}}_2 \hat{\boldsymbol{z}}_1^{{{\dagger}}} \rangle+ \langle \hat{\boldsymbol{z}}_1 \hat{\boldsymbol{z}}_2^{{{\dagger}}} \rangle + \langle \hat{\boldsymbol{z}}_2 \hat{\boldsymbol{z}}_2^{{{\dagger}}} \rangle , \end{gather}$$

(2.29)$$\begin{gather}\langle \hat{\boldsymbol{a}} \hat{\boldsymbol{a}} ^{{{\dagger}}} \rangle =\hat{\boldsymbol{\boldsymbol{\varGamma}}} \langle \hat{\boldsymbol{y}}_1 \hat{\boldsymbol{y}}_1^{{{\dagger}}} \rangle\hat{\boldsymbol{\boldsymbol{\varGamma}}}^{{{\dagger}}}. \end{gather}$$

Before further expanding these terms, we define

(2.30a,b)\begin{equation} \hat{{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{l}}}}}}(\omega) = {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{z}\boldsymbol{a}}^{{{\dagger}}}(\omega) {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{z}\boldsymbol{a}}(\omega) +{{\boldsymbol{\mathsf{P}}}} , \quad \hat{{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{r}}}}}}(\omega) ={-}{{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{z}\boldsymbol{a}}^{{{\dagger}}}(\omega) . \end{equation}

As shown in Appendix B, these terms are the counterparts of $\hat {{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{l}}}}}}$ and $\hat {{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{r}}}}}}$ for the full-knowledge control problem. It is now straightforward to show that

(2.31)$$\begin{gather} \langle \hat{\boldsymbol{y}}_1 \hat{\boldsymbol{y}}_1 ^{{{\dagger}}} \rangle = \hat{{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{l}}}}}}, \end{gather}$$

(2.32)$$\begin{gather}\langle \hat{\boldsymbol{z}}_1 \hat{\boldsymbol{z}}_1 ^{{{\dagger}}} \rangle = {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{z}\boldsymbol{f}} \hat{{{\boldsymbol{\mathsf{F}}}}}{{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{z}\boldsymbol{f}}^{{{\dagger}}}, \end{gather}$$

(2.33)$$\begin{gather}\langle \hat{\boldsymbol{z}}_2 \hat{\boldsymbol{z}}_2 ^{{{\dagger}}} \rangle = \hat{{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{r}}}}}}^{{{\dagger}}}\hat{\boldsymbol{\boldsymbol{\varGamma}}} \hat{{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{l}}}}}}\hat{\boldsymbol{\boldsymbol{\varGamma}}}^{{{\dagger}}} \hat{{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{r}}}}}}, \end{gather}$$

(2.34)$$\begin{gather}\langle \hat{\boldsymbol{z}}_2 \hat{\boldsymbol{z}}_1 ^{{{\dagger}}} \rangle ={-}\hat{{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{r}}}}}}^{{{\dagger}}} \hat{\boldsymbol{\boldsymbol{\varGamma}}} \hat{{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{r}}}}}}^{{{\dagger}}}. \end{gather}$$

The cost functional can then be expressed as

(2.35)\begin{align} J &= \int_{-\infty}^{\infty} {\rm Tr}( {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{z}\boldsymbol{f}} \hat{{{\boldsymbol{\mathsf{F}}}}}{{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{z}\boldsymbol{f}}^{{{\dagger}}} + \hat{{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{r}}}}}}^{{{\dagger}}}\hat{\boldsymbol{\boldsymbol{\varGamma}}}_{nc} \hat{{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{l}}}}}}\hat{\boldsymbol{\boldsymbol{\varGamma}}}_{nc}^{{{\dagger}}} \hat{{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{r}}}}}} - \hat{{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{r}}}}}}^{{{\dagger}}} \hat{\boldsymbol{\boldsymbol{\varGamma}}}_{nc} \hat{{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{r}}}}}}^{{{\dagger}}} {{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol{z}}}^{{{\dagger}}} - {{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol{z}}} \hat{{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{r}}}}}} \hat{\boldsymbol{\boldsymbol{\varGamma}}}_{nc}^{{{\dagger}}} \hat{{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{r}}}}}} )\, {\rm d} \omega \nonumber\\ &\quad + \int_{-\infty}^{\infty}{\rm Tr}( {{\boldsymbol{\mathsf{P}}}} \hat{\boldsymbol{\boldsymbol{\varGamma}}}_{nc} \hat{{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{l}}}}}} \hat{\boldsymbol{\boldsymbol{\varGamma}}}_{nc}^{{{\dagger}}} )\, {\rm d} \omega, \end{align}

where the subscript ‘nc’ was added to emphasize the non-causal nature of the control that will be obtained. Using the cyclic property of the trace and isolating terms with $\hat {\boldsymbol {\boldsymbol {\varGamma }}}_{nc}^{{\dagger} }$ , we obtain

(2.36)\begin{equation} J = \int_{-\infty}^{\infty} ( {\rm Tr} ( {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{z}\boldsymbol{f}} \hat{{{\boldsymbol{\mathsf{F}}}}}{{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{z}\boldsymbol{f}}^{{{\dagger}}} - \hat{{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{r}}}}}}^{{{\dagger}}} \hat{\boldsymbol{\boldsymbol{\varGamma}}}_{nc} \hat{{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{r}}}}}}^{{{\dagger}}} {{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol{z}}}^{{{\dagger}}} ) +{\rm Tr}(( \hat{{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{l}}}}}} \hat{\boldsymbol{\boldsymbol{\varGamma}}}_{nc} \hat{{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{l}}}}}} -\hat{{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{r}}}}}} {{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol{z}}} \hat{{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{r}}}}}} ) \hat{\boldsymbol{\boldsymbol{\varGamma}}}_{nc}^{{{\dagger}}} )) \, {\rm d} \omega, \end{equation}

where $\hat {{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{l}}}}}} = \hat {{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{r}}}}}}\hat {{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{r}}}}}}^{{\dagger} } + {{\boldsymbol{\mathsf{P}}}}$ was used. The minimum of $J$ is found by differentiating (2.36) with respect to $\hat {\boldsymbol {\boldsymbol {\varGamma }}}_{nc}^{{\dagger} }$, leading to

(2.37)\begin{equation} \hat{{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{l}}}}}} \hat{\boldsymbol{\boldsymbol{\varGamma}}}_{nc} \hat{{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{l}}}}}} = \hat{{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{r}}}}}} \hat{{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{r}}}}}}, \end{equation}

which can be solved for the optimal control kernel as

(2.38)\begin{equation} \hat{\boldsymbol{\boldsymbol{\varGamma}}}_{nc} = \hat{{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{l}}}}}}^{{-}1}\hat{{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{r}}}}}} \hat{{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{r}}}}}} \hat{{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{l}}}}}}^{{-}1} . \end{equation}

From (2.38), it can be seen that the optimal non-causal control kernel is a combination of the optimal non-causal estimation of targets ($\hat {{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{r}}}}}} \hat {{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{l}}}}}}^{-1}$, solution of (2.18)) with the optimal full-knowledge, non-causal control, described in Appendix B, ($\hat {{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{l}}}}}}^{-1}\hat {{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{r}}}}}}$, solution of (B3)), acting to minimize these targets with actuation inputs.

As this control kernel is, in general, not causal, it cannot be used in real-time applications. It does, however, provide upper bounds for the effectiveness of causal control, and it is the basis for the control method proposed by Sasaki et al. (Reference Sasaki, Tissot, Cavalieri, Silvestre, Jordan and Biau2018b), where the kernel is truncated to its causal part. This approach was also successfully applied to experiments (Brito et al. Reference Brito, Morra, Cavalieri, Araújo, Henningson and Hanifi2021; Maia et al. Reference Maia, Jordan, Cavalieri, Martini, Sasaki and Silvestre2021). However, truncating the non-causal control kernel is, in general, sub-optimal.

To obtain the optimal causal control law, causality is again enforced via Lagrange multipliers. The modified cost functional reads as

(2.39)\begin{align} J'& = J +\int_{-\infty}^{\infty} {\rm Tr}( \boldsymbol{\varLambda}_- (t) {\boldsymbol{\varGamma}}_c(t) + \boldsymbol{\varLambda}_-^{{{\dagger}}} (t) {\boldsymbol{\varGamma}}^{{{\dagger}}}_c(t) ) \, {\rm d} t \nonumber\\ & =J + \int_{-\infty}^{\infty} {\rm Tr}( \hat{\boldsymbol{\varLambda}}_- \hat{\boldsymbol{\boldsymbol{\varGamma}}}_c + \hat{\boldsymbol{\varLambda}}_-^{{{\dagger}}} \hat{\boldsymbol{\boldsymbol{\varGamma}}}^{{{\dagger}}}_c ) \,{\rm d} \omega, \end{align}

where the subscript $c$ emphasizes the causal nature of the control that will be obtained. Just as in the estimation problem, the linear terms added to $J$, contribute to an additional term when taking the derivative of $J'$. The control kernel is now given by the solution of

(2.40)\begin{equation} \hat{{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{l}}}}}} \hat{\boldsymbol{\boldsymbol{\varGamma}}}_c \hat{{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{l}}}}}} +\hat{\boldsymbol{\varLambda}}_- = \hat{{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{r}}}}}} {{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol{z}}} \hat{{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{r}}}}}}. \end{equation}

Again, as in the estimation problem, the requirements that ${\boldsymbol {\varGamma }(t<0)=0}$ and ${\boldsymbol {\varLambda }_-(t>0)=0}$ makes this a well-posed problem, and specifically another type of Wiener–Hopf problem. The procedure to solve this problem is equivalent to the one used for the estimation problem, and is presented in Appendix A.

The control kernel $\boldsymbol {\varGamma }'$ can be now recovered from $\boldsymbol {\varGamma }$. The expression for the actuation

(2.41)\begin{equation} \hat{\boldsymbol{a}} = \hat{\boldsymbol{\boldsymbol{\varGamma}}} \boldsymbol{y}_1 = \hat{\boldsymbol{\boldsymbol{\varGamma}}} (\hat{\boldsymbol{y}} -\hat{\boldsymbol{y}}_2) =\hat{\boldsymbol{\boldsymbol{\varGamma}}} (\hat{\boldsymbol{y}} -{{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}\boldsymbol{a}} \hat{\boldsymbol{a}}) , \end{equation}

can be rewritten as

(2.42)\begin{equation} \hat{\boldsymbol{a}} = \underbrace{ ({{\boldsymbol{\mathsf{I}}}} + \hat{\boldsymbol{\boldsymbol{\varGamma}}} {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}\boldsymbol{a}} )^{{-}1} \hat{\boldsymbol{\boldsymbol{\varGamma}}}}_{\hat{\boldsymbol{\boldsymbol{\varGamma}}}'} \hat{\boldsymbol{y}}, \end{equation}

thus recovering $\boldsymbol {\varGamma }'$.

The closed-loop control diagram is illustrated in figure 3, where the relation between $\hat {\boldsymbol {\boldsymbol {\varGamma }}}'$ and $\hat {\boldsymbol {\boldsymbol {\varGamma }}}$ is shown. As $\boldsymbol {y}_2$ is a function only of the previous actuation, it can be computed in real time and subtracted from $\boldsymbol {y}$ to obtain $\boldsymbol {y}_1$. This procedure is by construction included in ${\boldsymbol {\varGamma }}'$.

Figure 3.

Closed-loop control diagram.

Note that it is the process of removing the actuator's response from the sensor readings that can lead to instabilities if the feedback is not accurately modelled (Belson et al. Reference Belson, Semeraro, Rowley and Henningson2013). For amplifier flows, the actuators are typically located downstream of the sensors, the feedback tends to be small, and the optimal control is thus robust (Schmid & Sipp Reference Schmid and Sipp2016). This explains the successful use of optimal control in many studies (Semeraro et al. Reference Semeraro, Bagheri, Brandt and Henningson2011; Barbagallo et al. Reference Barbagallo, Dergham, Sipp, Schmid and Robinet2012; Morra et al. Reference Morra, Sasaki, Hanifi, Cavalieri and Henningson2020; Sasaki et al. Reference Sasaki, Morra, Cavalieri, Hanifi and Henningson2020; Maia et al. Reference Maia, Jordan, Cavalieri, Martini, Sasaki and Silvestre2021). As our method targets amplifier flows, the issue of robustness is not further addressed.

2.4. Recovering Kalman and LQR control gains

We now demonstrate that gain matrices for the Kalman-filter estimation (${{\boldsymbol{\mathsf{L}}}}$) and LQR control (${{\boldsymbol{\mathsf{K}}}}$) can be recovered from the Wiener–Hopf formalism. As both methods have the same optimality properties, they amount to different approaches for obtaining the same result if applied to a system satisfying the common assumptions in the derivation of Kalman filters and LQR control, such as white-in-time disturbances, i.e. $\hat {{{\boldsymbol{\mathsf{F}}}}}(\omega ) = {{\boldsymbol{\mathsf{F}}}}_0$.

From the Kalman-filter estimation equation (Åström & Wittenmark Reference Åström and Wittenmark2013),

(2.43)\begin{equation} \dfrac{{\rm d} \tilde{\boldsymbol{u}}}{{\rm d} t} = ({{\boldsymbol{\mathsf{A}}}} -{{\boldsymbol{\mathsf{L}}}} {{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol{y}}})\tilde{\boldsymbol{u}} + {{\boldsymbol{\mathsf{L}}}} \boldsymbol{y}, \end{equation}

the impulse response of the $i$th sensor, $y_j(t) =\delta _{j,i}\delta (t)$, is such that $\tilde {\boldsymbol {u}}(0^{+}) = {{\boldsymbol{\mathsf{L}}}}_i$, where ${{\boldsymbol{\mathsf{L}}}}_i$ is the $i$th row of ${{\boldsymbol{\mathsf{L}}}}$. Thus, from (2.16a,b),

(2.44)\begin{equation} {{\boldsymbol{\mathsf{L}}}} ={{\boldsymbol{\mathsf{T}}}}_{\boldsymbol{u}}(0). \end{equation}

The LQR control gains can be recovered by emulating an initial condition $\boldsymbol {u}_0$ at $t=0$, which can be accomplished by setting ${{{\boldsymbol{\mathsf{B}}}}_{\boldsymbol {f}}}={{\boldsymbol{\mathsf{I}}}}$ and using $\boldsymbol {f}(t) = \boldsymbol {u}_0 \delta (t)$, or equivalently, $\hat {\boldsymbol {z}}_1 = {{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol {z}}}{{\boldsymbol{\mathsf{R}}}} \boldsymbol {u}_0$. From the LQR framework,

(2.45)\begin{equation} \boldsymbol{a}(t) = {{\boldsymbol{\mathsf{K}}}}\boldsymbol{u}(t). \end{equation}

Since $\boldsymbol {u}(t<0)=0$, it follows that $\boldsymbol {a}(t<0) =0$. Comparing with the solution of (A16), (B6) and (2.10), we have

(2.46)\begin{align} \hat{\boldsymbol{a}}(\omega) & = \hat{{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{l}}}}}}_+^{{-}1}(\omega) (\hat{{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{l}}}}}}_-^{{-}1} (\omega) \hat{{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{r}}}}}} (\omega) \hat{\boldsymbol{z}}(\omega) )_+ \nonumber\\ & = \underbrace{\hat{{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{l}}}}}}_+^{{-}1}(\omega) (\hat{{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{l}}}}}}_-^{{-}1} (\omega) \hat{{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{r}}}}}} (\omega) {{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol{z}}} {{\boldsymbol{\mathsf{R}}}}(\omega) )_+}_{\hat{\boldsymbol{\varPi}}_c(\omega)} \boldsymbol{u}_0. \end{align}

Comparing (2.45) and (2.46), we have ${{\boldsymbol{\mathsf{K}}}} = \boldsymbol {\varPi }_c(0)$.

Knowledge of LQR gains ${{\boldsymbol{\mathsf{K}}}}$ and the Kalman filter ${{\boldsymbol{\mathsf{L}}}}$ may be useful to understand regions of the flow that require accurate estimation to obtain effective control strategies, as discussed by Freire et al. (Reference Freire, Cavalieri, Silvestre, Hanifi and Henningson2020).

2.5. Validation using a linearized Ginzburg–Landau problem

We here compare the kernels and gains computed with the present method to those obtained with the Kalman-filter and the LQR approaches. We use a linearized Ginzburg–Landau problem, for which standard tools can be used for the computation of estimation and control gains without resorting to model reduction. Such models have been used in many previous studies (Bagheri et al. Reference Bagheri, Henningson, Hœpffner and Schmid2009b; Lesshafft Reference Lesshafft2018; Cavalieri, Jordan & Lesshafft Reference Cavalieri, Jordan and Lesshafft2019). The model reads as

(2.47)\begin{equation} \frac{\partial \boldsymbol{u}(x,t)}{\partial t} = {{\boldsymbol{\mathsf{A}}}} \boldsymbol{u}(x,t) + \boldsymbol{f}(x,t) , \quad {{\boldsymbol{\mathsf{A}}}} ={-} U \frac{\partial }{\partial x} + \mu(x) + \gamma \frac{\partial^{2}}{\partial x^{2}},\end{equation}

and we use the same parameters as in Martini et al. (Reference Martini, Cavalieri, Jordan, Towne and Lesshafft2020), namely: $U=6 , \gamma =1-\mathrm {i}$ and $\mu (x)=\beta \mu _c(1-x/20)$, where $\beta =0.1$ was used and $\mu _c=U^{2} \textrm {Re} (\gamma ) / |\gamma |^{2}$ is the critical value for onset of absolute instability (Bagheri et al. Reference Bagheri, Henningson, Hœpffner and Schmid2009b).

For this validation, we use $\hat {{{\boldsymbol{\mathsf{F}}}}}={{\boldsymbol{\mathsf{I}}}}$, two sensors located at $x=5$ and $20$, an actuator at $x=15$, and a target at $x=30$. Figure 4 compares gains obtained from the proposed approach and from the algebraic Riccati equations. Both approaches produce identical gains, indicating their equivalence when white-noise forcing is assumed. Estimation and control kernels are shown in figure 5, again showing the equivalence between the two methods. The Kalman and LQG kernels are obtained as follows. From (2.16a,b), the state estimation is obtained as

(2.48)\begin{equation} \tilde {\boldsymbol{u}}(t) = \int_{-\infty}^{\infty} {{\boldsymbol{\mathsf{T}}}}_{u,kal}(\tau) \boldsymbol{y}(t-\tau) \, {\rm d} \tau, \end{equation}

with the estimation kernel given by

(2.49)\begin{equation} {{\boldsymbol{\mathsf{T}}}}_{u,kal}(\tau) = \begin{cases} {\rm e}^{ ({{\boldsymbol{\mathsf{A}}}}- {{\boldsymbol{\mathsf{L}}}} {{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol{y}}})\tau} {{\boldsymbol{\mathsf{L}}}}, & \tau\ge 0, \\ 0, & \tau<0. \end{cases} \end{equation}

The LQG kernel is obtained from

(2.50)$$\begin{gather} \frac{{\rm d}}{{\rm d} t} \tilde{\boldsymbol{u}}(t) = {{\boldsymbol{\mathsf{A}}}} \tilde {\boldsymbol{u}}(t) + {{{\boldsymbol{\mathsf{B}}}}_{\boldsymbol{a}}} \boldsymbol{a}(t)+{{\boldsymbol{\mathsf{L}}}} ( \boldsymbol{y}(t )- {{\boldsymbol{\mathsf{L}}}} {{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol{y}}} \tilde {\boldsymbol{u}}(t) ) , \end{gather}$$

(2.51)$$\begin{gather}\boldsymbol{a}(t) ={-} {{\boldsymbol{\mathsf{K}}}}\tilde {\boldsymbol{u}}(t), \end{gather}$$

with the solution reading as

(2.52)\begin{equation} \boldsymbol{a}(t) = \int_{-\infty}^{\infty} {\boldsymbol{\varGamma}}'_{lqg}(\tau) \boldsymbol{y}(t-\tau) \, {\rm d} \tau, \end{equation}

and where the control kernel is given by

(2.53)\begin{equation} {\boldsymbol{\varGamma}}'_{lqg}(\tau) = \begin{cases} {{\boldsymbol{\mathsf{K}}}} \, {\rm e}^{ ({{\boldsymbol{\mathsf{A}}}}-{{\boldsymbol{\mathsf{L}}}}{{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol{y}}}-{{{\boldsymbol{\mathsf{B}}}}_{\boldsymbol{a}}}{{\boldsymbol{\mathsf{K}}}}) \tau} {{\boldsymbol{\mathsf{L}}}}, & \tau\ge 0, \\ 0, & \tau<0. \end{cases} \end{equation}

Figures 4 and 5 show that the Riccati-based and proposed approaches provide the same results, illustrating their equivalence.

Figure 4.

(a) Blue/red lines show the estimation gains corresponding to the first and second sensor obtained using the proposed method. Dashed black lines show results obtained from the algebraic Riccati equation. (b) Same for control gains.

Figure 5.

(a,b) Kernels for the estimation of the system at different locations. Colour lines/black markers show results from the proposed method and from the algebraic Riccati equations. (c) Same for the control kernel. (d) Spatial support of the sensors, actuator and target. Vertical lines show locations for which the estimation kernels in (a,b) were computed.

3.1. Matrix-free implementation

In this section we apply methods developed in previous works (Martini et al. Reference Martini, Cavalieri, Jordan, Towne and Lesshafft2020, Reference Martini, Rodríguez, Towne and Cavalieri2021) to construct the terms $\hat {{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{l}}}}}} , \hat {{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{r}}}}}} , \hat {{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{l}}}}}}$ and $\hat {{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{r}}}}}}$, circumventing the need to construct the resolvent operator, which would make their construction prohibitive in any practical scenario. The approach consists of using time-domain solutions of the linearized equations to obtain the action of the resolvent operator on a vector, which in turn can be used to reconstruct the operator.

As an example, assume that the action of ${{\boldsymbol{\mathsf{R}}}}_{\boldsymbol {y}\boldsymbol {f}}$ on a vector $\hat {\boldsymbol {v}}_d$ can be efficiently computed. The operator can be constructed row-wise by setting $\hat {\boldsymbol {v}}_d=\boldsymbol {e}_i$, where $\boldsymbol {e}_i$ is an element of the canonical basis for the forcing space: the resulting vector ${{\boldsymbol{\mathsf{R}}}}_{\boldsymbol {y}\boldsymbol {f}}\boldsymbol {e}_i$ provides the $i$th row of ${{\boldsymbol{\mathsf{R}}}}_{\boldsymbol {y}\boldsymbol {f}}$. If $n_f$ is small, ${{\boldsymbol{\mathsf{R}}}}_{\boldsymbol {y}\boldsymbol {f}}$ can be recovered by repeating the procedure for $i=1,\ldots, n_f$, and subsequently used for the construction of $\hat {{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{l}}}}}}$.

To obtain the action of the resolvent operator on a vector for all resolved frequencies simultaneously, we use an approach based on the transient-response method, developed in Martini et al. (Reference Martini, Rodríguez, Towne and Cavalieri2021). The action of ${{\boldsymbol{\mathsf{R}}}}_{\boldsymbol {y}\boldsymbol {f}}$ and ${{\boldsymbol{\mathsf{R}}}}_{\boldsymbol {z}\boldsymbol {f}}$ on a vector $\hat {\boldsymbol {v}}_d(\omega )$ is obtained using the system

(3.1a–c)$$\begin{gather} \frac{{\rm d} \boldsymbol{q}}{{\rm d} t}(t) = {{\boldsymbol{\mathsf{A}}}} \boldsymbol{q}(t) + {{{\boldsymbol{\mathsf{B}}}}_{\boldsymbol{f}}} \boldsymbol{v}_d(t), \quad \boldsymbol{r}_d(t) = {{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol{y}}} \boldsymbol{q}(t),\quad \boldsymbol{s}_d(t) = {{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol{z}}} \boldsymbol{q}(t), \end{gather}$$

(3.2a–c)$$\begin{gather}\hat{\boldsymbol{q}} = {{\boldsymbol{\mathsf{R}}}}{{{\boldsymbol{\mathsf{B}}}}_{\boldsymbol{f}}} \hat{\boldsymbol{v}}_d, \quad \hat{\boldsymbol{r}}_d = \underbrace{{{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol{y}}}{{\boldsymbol{\mathsf{R}}}}{{{\boldsymbol{\mathsf{B}}}}_{\boldsymbol{f}}}}_{{{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}\boldsymbol{f}}} \hat{\boldsymbol{v}}_d, \quad \hat{\boldsymbol{s}}_d = \underbrace{{{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol{z}}}{{\boldsymbol{\mathsf{R}}}}{{{\boldsymbol{\mathsf{B}}}}_{\boldsymbol{f}}}}_{{{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{z}\boldsymbol{f}}} \hat{\boldsymbol{v}}_d, \end{gather}$$

where (3.2a–c) is obtained from a Fourier transform of (3.1a–c). Here $\hat {\boldsymbol {v}}_d$ is regarded as an input, with $\hat {\boldsymbol {r}}_d$ and $\hat {\boldsymbol {s}}_d$ as outputs of the system, with their dimensions implicit by the context. Using $\boldsymbol {v}_d(t)=\boldsymbol {e}_i\delta (t)$ ensures that the canonical basis is used for all frequencies. A time-domain solution of (3.1a–c) will be referred to as a forcing direct run. An actuator direct run is obtained by replacing ${{{\boldsymbol{\mathsf{B}}}}_{\boldsymbol {f}}}$ with ${{{\boldsymbol{\mathsf{B}}}}_{\boldsymbol {a}}}$, and provides the action of ${{\boldsymbol{\mathsf{R}}}}_{\boldsymbol {y}\boldsymbol {a}}$ and ${{\boldsymbol{\mathsf{R}}}}_{\boldsymbol {z}\boldsymbol {a}}$ on a given vector.

Similarly, the action of the operators ${{\boldsymbol{\mathsf{R}}}}_{\boldsymbol {y}\boldsymbol {f}}^{{\dagger} }$ and ${{\boldsymbol{\mathsf{R}}}}_{\boldsymbol {y}\boldsymbol {a}}^{{\dagger} }$ on a vector is constructed by time marching the adjoint system,

(3.3a–c)$$\begin{gather} -\frac{{\rm d} \boldsymbol{q}}{{\rm d} t}(t) = {{\boldsymbol{\mathsf{A}}}}^{{{\dagger}}} \boldsymbol{q}(t) + {{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol{y}}}^{{{\dagger}}} \boldsymbol{v}_a(t), \quad \boldsymbol{r}_a(t) = {{{\boldsymbol{\mathsf{B}}}}_{\boldsymbol{f}}}^{{{\dagger}}} \boldsymbol{q}(t), \quad \boldsymbol{s}_a(t) = {{{\boldsymbol{\mathsf{B}}}}_{\boldsymbol{a}}}^{{{\dagger}}} \boldsymbol{q}(t), \end{gather}$$

(3.4a–c)$$\begin{gather}\hat{\boldsymbol{q}} = {{\boldsymbol{\mathsf{R}}}}^{{{\dagger}}}{{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol{y}}}^{{{\dagger}}} \hat{\boldsymbol{v}}_a, \quad \hat{\boldsymbol{r}}_a = \underbrace{{{{\boldsymbol{\mathsf{B}}}}_{\boldsymbol{f}}}^{{{\dagger}}}{{\boldsymbol{\mathsf{R}}}}^{{{\dagger}}}{{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol{y}}}^{{{\dagger}}}}_{{{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}\boldsymbol{f}}^{{{\dagger}}}} \hat{\boldsymbol{v}}_a, \quad \hat{\boldsymbol{s}}_a = \underbrace{{{{\boldsymbol{\mathsf{B}}}}_{\boldsymbol{a}}}^{{{\dagger}}}{{\boldsymbol{\mathsf{R}}}}^{{{\dagger}}}{{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol{y}}}^{{{\dagger}}}}_{{{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}\boldsymbol{a}}^{{{\dagger}}}} \hat{\boldsymbol{v}}_a, \end{gather}$$

referred to here as a sensor adjoint run. A target adjoint run is obtained by replacing ${{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol {y}}}$ with ${{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol {z}}}$ and can be used to obtain the action of ${{\boldsymbol{\mathsf{R}}}}_{\boldsymbol {z}\boldsymbol {f}}^{{\dagger} }$ and ${{\boldsymbol{\mathsf{R}}}}_{\boldsymbol {z}\boldsymbol {a}}^{{\dagger} }$.

The action of more complex terms can be obtained by solving (3.1a–c) and (3.3a–c) in succession. For example, if the output $\boldsymbol {r}_a$ of a sensor adjoint run is an input of a direct run, then the outputs of the latter will be given by ${ \hat {\boldsymbol {r}}_d = {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol {y}\boldsymbol {f}}{{\boldsymbol{\mathsf{R}}}}_{\boldsymbol {y}\boldsymbol {f}}^{{\dagger} }\hat {\boldsymbol {v}}_a}$ and ${ \hat {\boldsymbol {s}}_d = {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol {z}\boldsymbol {f}}{{\boldsymbol{\mathsf{R}}}}_{\boldsymbol {y}\boldsymbol {f}}^{{\dagger} }\hat {\boldsymbol {v}}_a }$. The resulting system is referred to as a sensor adjoint-direct run. Following similar procedures, direct-adjoint, direct-adjoint-direct and adjoint-direct-adjoint runs are constructed. The terms whose action are obtained from each of these runs are illustrated in figure 6.

Figure 6.

Illustration of the terms obtained from different combinations of direct and adjoint time-domain solutions. The left and right flowcharts represent a sensor adjoint(-direct-adjoint) system and an actuator direct(-adjoint-direct) system. Readings refer to inner products of the state with sensors/forcing/targets/actuators spatial supports, and F.T. to a Fourier transform. Diagrams for the target and forcing systems are equivalent to the left(right) diagrams with ${{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol {y}}}({{{\boldsymbol{\mathsf{B}}}}_{\boldsymbol {a}}})$ exchanged by ${{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol {z}}}({{{\boldsymbol{\mathsf{B}}}}_{\boldsymbol {f}}})$, and the $\boldsymbol {y}(\textit {a})$ subscripts of the rightmost terms in each quantity exchanged by $\boldsymbol {z}(f)$.

In practice, to link together direct and adjoint runs as described above, checkpoints of the output of one time integration are saved to disk and subsequently read and interpolated in the following run to be used as forcing terms. Details of this procedure as well as the impact of different interpolation strategies are described by Martini et al. (Reference Martini, Rodríguez, Towne and Cavalieri2021), to which we refer the reader for details. Alternatively, Farghadan et al. (Reference Farghadan, Towne, Martini and Cavalieri2021) developed an approach that minimizes the data to be retained using streaming Fourier sums to obtain the action of the operator on a discrete set of frequencies.

The most effective approach for constructing the Wiener–Hopf problems depend on the rank of forcing and targets. In the following sections, we outline the best approach for four possible scenarios and discuss the associated computational cost and the types of parametric studies that can be easily conducted in each case.

3.1.5. Summary

A brief summary of the costs of each scenario is presented in table 1. Note that the descriptions above focused on the solution of the Wiener–Hopf problem, which provide $\boldsymbol {\varGamma }$. In some scenarios, supplementary runs are required to obtain ${{\boldsymbol{\mathsf{R}}}}_{\boldsymbol {y}\boldsymbol {a}}$, which is required to construct $\boldsymbol {\varGamma }'$ as in (2.42). This term can be obtained using sensor adjoint runs or actuator direct runs. We also report the extra runs required to obtain ${{\boldsymbol{\mathsf{R}}}}_{\boldsymbol {z}\boldsymbol {f}}^{{\dagger} }{{\boldsymbol{\mathsf{R}}}}_{\boldsymbol {z}\boldsymbol {f}}$, which will be used to obtain offline estimates of the control performance in § 3.5.

Table 1.

Description and number of time-domain solutions needed in each scenario to construct the $\boldsymbol {\varGamma }$, and the extra time-domain solutions required to construct $\boldsymbol {\varGamma }'$ and to perform an offline estimation of the control performance.

3.2. Using coloured forcing statistics

Previous results show that incorporating accurate coloured forcing statistics in the model is important to obtain accurate estimates of the flow state (Chevalier et al. Reference Chevalier, Hœpffner, Bewley and Henningson2006; Martini et al. Reference Martini, Cavalieri, Jordan, Towne and Lesshafft2020; Amaral et al. Reference Amaral, Cavalieri, Martini, Jordan and Towne2021). We now detail how to incorporate forcing colour for each of the scenarios described in the previous sections.

In the scenarios described in §§ 3.1.1 and 3.1.3, since the term ${{\boldsymbol{\mathsf{R}}}}_{\boldsymbol {y}\boldsymbol {f}}$ is a small matrix, forcing colour could be easily included a posteriori in the construction of $\hat {{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{l}}}}}}$ and $\hat {{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{r}}}}}}$. In the scenarios §§ 3.1.2 and 3.1.4, where ${{\boldsymbol{\mathsf{R}}}}_{\boldsymbol {y}\boldsymbol {f}}{{\boldsymbol{\mathsf{R}}}}_{\boldsymbol {y}\boldsymbol {f}}^{{\dagger} }$ is directly obtained from a time-domain solution, the inclusion of forcing colour requires time-domain convolutions of the output of the sensor adjoint system with the forcing cross-correlation matrices, before its use as an input of the direct system. However, the convolutions of a large matrix (${{\boldsymbol{\mathsf{F}}}}$) and a vector (output of the adjoint system) is typically unfeasible. To circumvent this limitation, the approach proposed by Martini et al. (Reference Martini, Cavalieri, Jordan, Towne and Lesshafft2020) can be used. An estimated forcing CSD, obtained from low-rank flow measurements, is used to construct $\hat {{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{l}}}}}}$ and $\hat {{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{r}}}}}}$. This can be done without explicit construction of $\hat {{{\boldsymbol{\mathsf{F}}}}}$, as will be shown next.

From a set of $n_{y'}$ auxiliary sensor readings, $\boldsymbol {y}'$, obtained as ${{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol {y}}}'\boldsymbol {u}$, which we assume to contain the set of sensors that will be used for estimation and/or control ($\boldsymbol {y}$), and their CSD, ${{\boldsymbol{\mathsf{Y}}}}'$, the estimated forcing colour is obtained as

(3.5)\begin{equation} \hat{{{\boldsymbol{\mathsf{F}}}}}' = \hat{\boldsymbol{\mathsf{T}}}_f \hat{\boldsymbol{\mathsf{Y}}}' \hat{\boldsymbol{\mathsf{T}}}_f^{\dagger}, \end{equation}

where

(3.6)\begin{equation} \hat{{{\boldsymbol{\mathsf{T}}}}}_f' ={{\boldsymbol{\mathsf{R}}}}^{{{\dagger}}}_{y'f} ( {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}'\boldsymbol{f}} {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}'\boldsymbol{f}}^{{{\dagger}}} + \hat{\boldsymbol{\mathsf{N}}} ')^{{-}1}, \end{equation}

is a transfer function that estimates forcing from measurements with a prior assumption of white-noise forcing (Martini et al. Reference Martini, Cavalieri, Jordan, Towne and Lesshafft2020; Towne et al. Reference Towne, Lozano-Durán and Yang2020).

Estimates for $\hat {{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{l}}}}}}$ and $\hat {{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{r}}}}}}$ using the estimated forcing CSD, $\hat {{{\boldsymbol{\mathsf{F}}}}}'$, referred to here as $\hat {{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{l}}}}}}'$ and $\hat {{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{r}}}}}} '$, read as

(3.7)$$\begin{gather} \hat{{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{l}}}}}}' = {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}\boldsymbol{f}} \hat{{{\boldsymbol{\mathsf{F}}}}}' {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}\boldsymbol{f}}^{{{\dagger}}} + \hat{\boldsymbol{\mathsf{N}}} = {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}\boldsymbol{f}} \hat{{{\boldsymbol{\mathsf{T}}}}}_f \hat{\boldsymbol{\mathsf{Y}}}' \hat{{{\boldsymbol{\mathsf{T}}}}}_f^{{{\dagger}}} {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}\boldsymbol{f}}^{{{\dagger}}} + \hat{\boldsymbol{\mathsf{N}}}, \end{gather}$$

(3.8)$$\begin{gather}\hat{{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{r}}}}}} ' = {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{z}\boldsymbol{f}} \hat{{{\boldsymbol{\mathsf{F}}}}}' {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}\boldsymbol{f}}^{{{\dagger}}} = {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{z}\boldsymbol{f}} \hat{{{\boldsymbol{\mathsf{T}}}}}_f \hat{\boldsymbol{\mathsf{Y}}}' \hat{{{\boldsymbol{\mathsf{T}}}}}_f^{{{\dagger}}} {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}\boldsymbol{f}}^{{{\dagger}}}, \end{gather}$$

Note that $\hat {{{\boldsymbol{\mathsf{T}}}}}' _f \in \mathbb {C}^{n_f \times n_{y'}}$, and that $n_f$ was assumed to be large. To avoid operations with large matrices, the compound term

(3.9)\begin{equation} {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{z}\boldsymbol{f}} \hat{{{\boldsymbol{\mathsf{T}}}}}' _f= {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{z}\boldsymbol{f}} {{\boldsymbol{\mathsf{R}}}}^{{{\dagger}}}_{\boldsymbol{y}'\boldsymbol{f}} ( {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}'\boldsymbol{f}} {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}'\boldsymbol{f}}^{{{\dagger}}} + \hat{\boldsymbol{\mathsf{N}}} ')^{{-}1}, \end{equation}

that has size $n_y \times n_{y'}$, can be constructed from ${ {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol {z}\boldsymbol {f}} {{\boldsymbol{\mathsf{R}}}}^{{\dagger} }_{\boldsymbol {y}'\boldsymbol {f}} \in \mathbb {C}^{n_z\times n_{y'}}}$ and ${{\boldsymbol{\mathsf{R}}}}_{\boldsymbol {y}'\boldsymbol {f}} {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol {y}'\boldsymbol {f}}^{{\dagger} } \in \mathbb {C}^{n_{y'}\times n_{y'}}$. As the terms extra terms ${{\boldsymbol{\mathsf{R}}}}_{\boldsymbol {y}'\boldsymbol {f}}$ and ${{\boldsymbol{\mathsf{R}}}}_{\boldsymbol {y}'\boldsymbol {f}} {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol {y}'\boldsymbol {f}}^{{\dagger} }$ can be obtained from the matrix-free approach described, the time-domain convolution can be replaced by a few extra time-domain solutions of the linearized problems and inexpensive post-processing.

3.3. Adjointless method

We now show that the terms needed to construct the estimation and control kernels can be obtained from experimental data alone. We first write the sensors and targets in terms of the external disturbances,

(3.10)\begin{equation} \left[\begin{matrix} \hat {\boldsymbol{y}}_1\\ \hat {\boldsymbol{z}}_1 \end{matrix}\right] = \left[\begin{matrix} {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}\boldsymbol{f}}\\\ {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{z}\boldsymbol{f}} \end{matrix}\right] \hat {\boldsymbol{f}} + \left[\begin{matrix} 1 \\0 \end{matrix}\right] \hat {\boldsymbol{n}}. \end{equation}

Their CSDs are then obtained from the forcing statistics (Towne et al. Reference Towne, Lozano-Durán and Yang2020) as

(3.11)\begin{equation} \left[\begin{matrix} {{{\boldsymbol{\mathsf{S}}}}}_{y_1,y_1} & {{{\boldsymbol{\mathsf{S}}}}}_{y_1,z_1} \\ {{{\boldsymbol{\mathsf{S}}}}}_{z_1,y_1} & {{{\boldsymbol{\mathsf{S}}}}}_{z_1,z_1} \end{matrix}\right] = \left[\begin{matrix} {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}\boldsymbol{f}}\\\ {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{z}\boldsymbol{f}} \end{matrix}\right] \hat {{\boldsymbol{\mathsf{F}}}} \left[\begin{matrix} {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}\boldsymbol{f}} ^{{{\dagger}}} & {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{z}\boldsymbol{f}}^{{{\dagger}}} \end{matrix}\right] + \left[\begin{matrix} \hat {{\boldsymbol{\mathsf{N}}}} & 0 \\ 0 & 0 \end{matrix}\right] =\left[\begin{matrix} \hat{{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{l}}}}}} & \hat{{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{r}}}}}} \\ \hat{{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{r}}}}}} ^{{{\dagger}}} & {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{z}\boldsymbol{f}} \hat {{\boldsymbol{\mathsf{F}}}} {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{z}\boldsymbol{f}}^{{{\dagger}}} \end{matrix}\right]. \end{equation}

The non-causal estimation is obtained via the transfer function ${{\boldsymbol{\mathsf{S}}}}_{y,z} {{\boldsymbol{\mathsf{S}}}}_{y,y}^{-1}$. The terms $\hat {{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{l}}}}}}$ and ${{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol {z}}} \hat {{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{r}}}}}}$ can be constructed from measurements of the uncontrolled system, as they correspond to sensor CSDs (${{{\boldsymbol{\mathsf{S}}}}}_{y_1,y_1}$) and the cross-spectra between sensors and targets (${{{\boldsymbol{\mathsf{S}}}}}_{y_1,z_1}$), respectively. As the terms $\hat {{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{l}}}}}}$ and $\hat {{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{r}}}}}}$ can be obtained from the actuators’ impulse response, all the required terms can be obtained using physical or numerical experiments. This not only provides a simple framework to obtain optimal control based on a data-driven approach, but it also relaxes the requirement of an adjoint solver to obtain these control strategies for large systems, either by using readings from the nonlinear problem, as done by Martinelli (Reference Martinelli2009), or from the direct linearized problem excited by stochastic forcing.

The adjointless approach presented here is in fact the classical application of the Wiener regulator, which is constructed directly from sensor/target CSDs, with the derivation presented here showing that the two approaches, the classical and the resolvent-based, are equivalent when the exact force model is used. However, the statistical convergence of the CSDs is considerably more expensive, and typically less accurate, than its construction from first principles (§ 3.1) and, thus, the latter is preferable to a numerical experiment if an adjoint solver is available. For physical experiments, where obtaining a long time series is typically inexpensive, the method presented here provides an efficient way to obtain a data-driven optimal control law.

3.4. The role of sensors for the estimation problem

Insights into the problem of sensor placement can be obtained from further analysis of (3.11), focusing the discussion on non-causal estimation and white-noise forcing for simplicity. Sensor and target sensitivities to the forcing terms are given by ${{\boldsymbol{\mathsf{R}}}}_{\boldsymbol {y}\boldsymbol {f}}^{{\dagger} }$ and ${{\boldsymbol{\mathsf{R}}}}_{\boldsymbol {z}\boldsymbol {f}}^{{\dagger} }$. We define two forcing subspaces, $\hat {{{\boldsymbol{\mathsf{W}}}}}_y$ and $\hat {{{\boldsymbol{\mathsf{W}}}}}_z$, spanned by the rows of ${{\boldsymbol{\mathsf{R}}}}_{\boldsymbol {y}\boldsymbol {f}}^{{\dagger} }$ and ${{\boldsymbol{\mathsf{R}}}}_{\boldsymbol {z}\boldsymbol {f}}^{{\dagger} }$. Defining $\hat {\boldsymbol {w}}_{i,y}$ and $\hat {\boldsymbol {w}}_{i,z}$ to be orthogonal bases for these subspaces, the forcing CSD can be decomposed into four different subspaces:

• • $\mathcal {F}_{y,z}$: spanned by components given by $\hat {\boldsymbol {w}}_{i,y}\hat {\boldsymbol {w}}_{i,z}^{{\dagger} } + \hat {\boldsymbol {w}}_{i,z}\hat {\boldsymbol {w}}_{i,y}^{{\dagger} }$, corresponds to forcing components correlated with responses in both $\boldsymbol {y}$ and $\boldsymbol {z}$,

• • $\mathcal {F}_{y}$: the subspace spanned by components $\hat {\boldsymbol {w}}_{i,y}\hat {\boldsymbol {w}}_{i,y}^{{\dagger} }$ that are orthogonal to $\mathcal {F}_{y,z}$, containing forcing components correlated with responses in $\boldsymbol {y}$ but not in $\boldsymbol {z}$,

• • $\mathcal {F}_{z}$: the subspace spanned by components $\hat {\boldsymbol {w}}_{i,z}\hat {\boldsymbol {w}}_{i,z}^{{\dagger} }$ that are orthogonal to $\mathcal {F}_{y,z}$, containing forcing components correlated with responses in $\boldsymbol {z}$ but not in $\boldsymbol {y}$,

• • $\mathcal {F}^{\perp }$: the complement of the three above subspaces, containing forcing components that are correlated with neither $\boldsymbol {y}$ nor $\boldsymbol {z}$.

Defining $\hat {{{\boldsymbol{\mathsf{F}}}}}_{y}$ as the projection of $\hat {{{\boldsymbol{\mathsf{F}}}}}$ into $\mathcal {F}_y$, with similar definitions for the other subspaces, (3.11) can be rewritten as

(3.12)\begin{align} \left[\begin{matrix} {{{\boldsymbol{\mathsf{S}}}}}_{y_1,y_1} & {{{\boldsymbol{\mathsf{S}}}}}_{y_1,z_1} \\ {{{\boldsymbol{\mathsf{S}}}}}_{z_1,y_1} & {{{\boldsymbol{\mathsf{S}}}}}_{z_1,z_1} \end{matrix}\right] &= \left[\begin{matrix} {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}\boldsymbol{f}} \hat {{\boldsymbol{\mathsf{F}}}}_y {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}\boldsymbol{f}}^{{{\dagger}}} +\hat {{\boldsymbol{\mathsf{N}}}} & 0\\ 0 & 0 \end{matrix}\right] \nonumber\\ &\quad +\left[\begin{matrix} {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}\boldsymbol{f}} \hat {{\boldsymbol{\mathsf{F}}}}_{y,z} {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}\boldsymbol{f}}^{{{\dagger}}} & {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}\boldsymbol{f}} \hat {{\boldsymbol{\mathsf{F}}}}_{y,z} {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{z}\boldsymbol{f}}^{{{\dagger}}}\\ {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{z}\boldsymbol{f}} \hat {{\boldsymbol{\mathsf{F}}}}_{y,z} {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{y}\boldsymbol{f}}^{{{\dagger}}} & {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{z}\boldsymbol{f}} \hat {{\boldsymbol{\mathsf{F}}}}_{y,z} {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{z}\boldsymbol{f}}^{{{\dagger}}} \end{matrix}\right] + \left[\begin{matrix} 0 & 0 \\ 0 & {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{z}\boldsymbol{f}} \hat {{\boldsymbol{\mathsf{F}}}}_z {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{z}\boldsymbol{f}}^{{{\dagger}}} \end{matrix}\right]. \end{align}

Each of these subspaces plays a different role in the estimation. The subspace $\mathcal {F}_y$ generates responses at the sensors but not at the target and, thus, it does not provide any useful information that can be used to construct the estimates, instead appearing in the problem in a similar way as does the sensor noise. On the other hand, the subspace $\mathcal {F}_{y,z}$ generates responses at sensors and targets and, thus, the response to this forcing measured by the sensors can be used to estimate the corresponding response of the flow at the target locations. This is the term that is effectively used for target estimation. Finally, the subspace $\mathcal {F}_{z}$ corresponds to forces that have responses at the target but not at the sensors. Accordingly, the responses associated with it cannot be estimated.

This analysis has a direct relation with the concept of observable forces discussed in previous works by Towne et al. (Reference Towne, Lozano-Durán and Yang2020): only target components that are excited by forcing components that also generate readings on the sensors can be estimated and, as a consequence, controlled. Martini et al. (Reference Martini, Cavalieri, Jordan, Towne and Lesshafft2020) discussed how the correlation between difference forcing components can allow the estimation of the responses to non-observable forcing components, and this was shown to considerably improve turbulent flow estimation. Note that, in general, non-causal estimation/control is needed for use of the full correlation, i.e. estimation/control of all the target components which are correlated with the sensors. The causality constraint allows for a real-time application at the price of deteriorating the estimation/control.

To optimize the non-causal estimation, one should seek to minimize ${{\boldsymbol{\mathsf{F}}}}_z$ and reduce the effect of ${{\boldsymbol{\mathsf{F}}}}_y$ on the estimation of ${{\boldsymbol{\mathsf{F}}}}_{y,z}$. Although causality imposes extra restrictions, which further complicates the problem, the above discussion can provide insights into the sensor placement problem, which is still an active topic of research.

4.1. Control of the linearized problem

We initially focus on a linearized system disturbed by high-rank external forces, designing control strategies to minimize readings from low-rank targets. We thus perform a limited study on the sensor placement, but a large actuator placement study. These results will be used to inform sensors and actuators to be used on a full-rank target scenario.

4.1.1. Control with a single sensor

We first consider the use of a single sensor (middle circle in figure 7) and investigate the placement of one actuator for control (locations indicated by crosses in figure 7). Three control approaches are explored: non-causal control (nc), the truncated non-causal control (tnc) and optimal-causal control (c). Truncated non-causal control corresponds to a truncation of the optimal non-causal control to its causal part, corresponding to the approach used in previous works (Sasaki et al. Reference Sasaki, Tissot, Cavalieri, Silvestre, Jordan and Biau2018b; Brito et al. Reference Brito, Morra, Cavalieri, Araújo, Henningson and Hanifi2021), which was seen, in some cases, to be a good approximation for optimal causal control, pointing to a wave-cancelling nature of optimal control strategies in these scenarios (Sasaki et al. Reference Sasaki, Morra, Fabbiane, Cavalieri, Hanifi and Henningson2018a). Non-causal control cannot be used in real-time applications, but it provides upper bounds for the performance of any linear control strategy and, thus, can be useful to evaluate different sensor and actuator placements.

As we use high-rank forcing and a low-rank target, the configuration corresponds to scenario (ii) described in § 3.1. A total of three time-domain solutions were obtained, corresponding to one adjoint-direct system for the sensors, from which the estimation terms ${{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{l}}}}}$ and ${{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{r}}}}}$ are obtained, and one adjoint-direct system for the targets, from which impulse responses from any actuator can be obtained. All the following results were obtained via inexpensive post-processing of the resulting data.

Representative control kernels for different actuator placements are shown in figure 8. The causal kernel is seen to converge to the non-causal kernel for large $\tau$ when the actuator is not far upstream of the sensor. When the non-causal control law uses non-causal information, the causal control law typically exhibits a spike at $\tau =0$. This behaviour can be understood as the control law using the most recent information to compensate for future information, which is not available. Similar behaviour was reported by Morra et al. (Reference Morra, Sasaki, Hanifi, Cavalieri and Henningson2020).

Figure 8.

Causal and non-causal control for different actuator placements and a sensor at $(x,y)=(-1,0)$. (a) Actuator at $x=-3.5$. (b) Actuator at $x= 0.5$. (c) Actuator at $x= 1.5$.

To compare the different control strategies, (3.13) is used to compute the target's power spectral density (PSD) and expected energy, i.e. the integral of the PSD, with the results presented in figure 9. All control strategies converge approximately to the same PSD reduction when the actuator is located downstream of the sensor. Not surprisingly, non-causal control provides the highest PSD reduction for all configurations, with little dependence on the actuator placement. As the non-causal control law can be solved independently for each frequency, the PSDs obtained with this control are lower than that of the uncontrolled case for all frequencies. Different trends are observed for the causal or the truncated-non-causal (TNC) control. Comparing the target PSDs for the causal and non-causal cases in figure 9(b), the non-causal control provides the lowest PSD at all frequencies. The causal control increases target PSD at lower frequencies ($|\omega |<0.2$), which is however compensated by the reduction for the other frequencies. Such a compromise between PSD reduction in different frequencies is unavoidable in causal control strategies, due to constraints imposed by causality. Figure 10 compares the kernels in the frequency domain and shows that each approach better represents the non-causal kernel for different frequencies. A comparison of figures 9(b) and 10 shows that the causal control better reproduces the non-causal control for the most relevant frequencies, i.e. those which provide a higher reduction in the target PSD, explaining its superior performance.

Figure 9.

Target reading PSDs and energy for different control strategies and actuator placement. In (a)–(c) and for the solid lines in (d), ${{\boldsymbol{\mathsf{P}}}}=10^{-2}\max (|{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{l}}}}}(\omega )|)$ was used. A value of ${{\boldsymbol{\mathsf{P}}}}=10^{-5}\max (|{{\boldsymbol{\mathsf{H}}}}_{{{\boldsymbol{\mathsf{l}}}}}(\omega )|)$ was used for dotted lines in (d). (a) Actuator at $x=-9.5$. (b) Actuator at $x=-3.5$. (c) Actuator at $x= 0.5$. (d) Target energy.

Figure 10.

Difference between causal, and TNC control kernels with respect to the non-causal kernel in the frequency domain. The configuration is the same as in figure 9(b). As the quantities are complex, the absolute value of the difference is shown to measure how the causal and TNC approaches deviate from the non-causal kernel.

Figure 9(d) shows that the three approaches show different trends as the actuator is moved upstream. While actuator placement has a negligible influence on the non-causal control, it significantly affects the TNC control, leading to an increase in the target energy for some configurations. Causal control, on the other hand, remains efficient even when located upstream of the sensor, with the reduction in the PSD degraded by $50\,\%$ only when the actuator is located approximately 4.5 step heights upstream of the sensor.

Two different effects can explain these observations. Structures that emerge from the upstream disturbances can have significant coherence lengths, of the order of two convective time units, as later reported in figure 15. This coherence length allows the sensor to partially estimate their upstream components, which can thus be cancelled by the actuator. It is speculated that a similar mechanism is responsible for the control of a jet (a convectively unstable flow, as the present example) using downstream sensors obtained by Bychkov et al. (Reference Faranosov, Belyaev, Bychkov, Kopiev, Kopiev, Moralev and Kazansky2019) and Kopiev et al. (Reference Kopiev, Faranosov, Kopiev, Bychkov, Moralev and Kazansky2019). An extreme example is when the flow is excited with a harmonic, rank-1, forcing: forces and flow response have infinite coherence lengths, and can thus be estimated and controlled based on only a few sensor readings, from which amplitudes and phases are extracted.

Another possible mechanism is the exploration of different control approaches by the causal control. As the flow studied here is incompressible, the actuation has an instantaneous, although possibly small, effect throughout the flow. This response is typically negligible for the non-causal control, which favours more effective mechanisms, but may be the only one available for the causal control when the actuator is located far upstream and exploited if low actuation penalties are used. This interpretation is supported by the dotted lines in figure 9, where the actuation penalty was drastically reduced. For upstream actuators, a significant reduction of the targets PSDs is observed for the causal control, but a small effect is observed for the non-causal control. Actuator placements upstream of the sensors tend to be less robust to unmodelled dynamics (Belson et al. Reference Belson, Semeraro, Rowley and Henningson2013), and it is a configuration that is typically not used for amplifier flows (Schmid & Sipp Reference Schmid and Sipp2016). Thus, upstream actuators are not further investigated in this work.

4.1.2. Control with multiple sensors

Hervé et al. (Reference Hervé, Sipp, Schmid and Samuelides2012) reported higher perturbation reductions for a similar control configuration than those obtained here using a single sensor and actuator. However, the cited work considered only rank-1 disturbances. The present work deals with a more complex forcing scenario, with several waves exciting the Kelvin–Helmholtz instability. The identification of such multiple waves invariably requires multiple sensors; we thus explore this scenario. Four additional time-domain solutions were necessary for the results presented next, corresponding to adjoint-direct systems for the added sensors, the top and bottom circles in figure 7.

Figure 11 shows the expected target energy using the three sensors. A comparison with results obtained using one sensor identifies two trends: (i) lower target PSDs are observed, indicating that the additional sensors are indeed necessary to identify the multiple incoming waves and to accurately estimate the target readings; (ii) a larger distance between the actuator and the sensors is required for TNC to reproduce the optimal control strategies. Comparing the kernels when using one and three sensors in figure 12, it can be seen that non-causal control for three sensors is more ‘spread’ in $\tau$ and, thus, has more signal content for $\tau <0$. It is then expected that its truncation leads to a more significant degradation of the control.

Figure 11.

Same as figure 9, but using three sensors. Results for SISO control are reproduced in the dashed lines for reference.

Figure 12.

Control kernels using an actuator at $x=-1.5$ and (a) one sensor, located at $y=0.5$, (b) the three sensors indicated in figure 7. (a) Control kernels using one sensor and (b) control kernels using three sensors.

4.1.3. Full-rank target control

From figures 9 and 11, we conclude that the actuator should be located shortly downstream of the sensors. With sensor and actuator placements defined, control kernels for full-rank targets are constructed. This required that the following additional time-domain solutions be performed: the last equation in the sensor direct-adjoint-direct problem for each of the sensors, and the full direct-adjoint system for each actuator. In total nine additional time-domain solutions were used.

Figure 13 compares the control kernel for a rank-1 and full-rank ${{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol {z}}}$. The close match between the two kernels indicates that downstream of the step the problem can be well approximated by a rank-1 model, and that the high-rank behaviour of the problem is indeed restricted to the receptivity of the Kelvin–Helmholtz modes to channel disturbances. This scenario is further supported by results presented in figure 14, where flow snapshots and time series for the perturbation norm are shown for several configurations of the sensor, actuators and targets. Control designed for full-rank targets and multiple actuators do lead to perturbation energy reduction, but the biggest improvement is provided by the use of multiple sensors.

Figure 13.

Comparison between control kernels for low- and full-rank targets.

Figure 14.

Snapshots for the uncontrolled flow (a) and controlled flows using the causal control kernel (b)–(g). The sensors and actuators used in each scenario are indicated by circles and crosses. Low-rank targets are indicated by triangles on the figures, full-rank targets, i.e. ${{{\boldsymbol{\mathsf{C}}}}_{\boldsymbol {z}}}= {{\boldsymbol{\mathsf{I}}}}$, were used when the figures do not contain triangles. In (h), norms for perturbations for $x>0$, with colours highlighting sensor and actuators used and solid/dotted lines low/full target ranks.

4.2. Adjoint-less and empirical application

Here we illustrate a data-driven application of the proposed method. The necessary data can be obtained from experimental set-ups, from numerical solutions of the nonlinear problem or from the direct linearized problem disturbed by stochastic forcing. We focus here on the latter scenario, with the use of a single sensor and target. The generalization to more complex configurations, with more sensors, actuators and targets, is straightforward. From sensor and target time series, the auto-correlation and cross-correlations are computed in the time domain. As previously described in § 2, (3.11), the terms $\hat {{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{l}}}}}}$ and $\hat {{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{r}}}}}}$ are equivalent to the corresponding CSDs.

Figure 15 compares estimates of ${{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{l}}}}}$ and ${{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{r}}}}}$, obtained ‘empirically’, from the numerical experiment with different data lengths, and ‘analytically’, obtained with time marching of direct and adjoint equations. The resulting control kernels are also compared. Convergence trends are shown using the error, defined as

(4.1)\begin{equation} error = \frac{\| \tilde{{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{l}}}}} }-{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{l}}}}}\| }{\|{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{l}}}}}\|}, \end{equation}

where the $L^{2}$ norm is used. The terms $\tilde {{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{l}}}}}}$ and ${{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{l}}}}} }$ indicate, respectively, the empirical and analytical terms, with similar expressions for the other terms. The trends suggest convergence rates scaling with $\Delta T^{-1/2}$, where $\Delta T$ is the time-series length used to estimate ${{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{l}}}}}$ and ${{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{r}}}}}$.

Figure 15.

Comparison between ${{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{l}}}}} , {{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{r}}}}} , {\boldsymbol {\varGamma }}_c$ and errors between the ‘empirical’ functions and those obtained using the time-march method proposed. Results are presented for several time lengths used when estimating the auto- and cross-correlations. Results with the time-march method require time integration of approximately 300 time units, indicated by the dotted line in (d). A white sensor noise was added to $\hat {{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{l}}}}}}$ to regularize its factorization.

Although significantly more expensive, and typically less accurate than results obtained with the use of adjoint solvers, the use of numerical experiments does extend the applicability of the method to virtually any scenario and solver, and also to derive control laws directly from experiments, which typically can be carried out to obtain long time series for better convergence of the required correlations. Note that if experimental data are used, sensor noise is already present in the data. However, if noise levels are small, adding a small noise may be required to better condition the factorization of $\hat {{{\boldsymbol{\mathsf{G}}}}_{{{\boldsymbol{\mathsf{l}}}}}}$.

4.3. Estimation and control of the nonlinear system

Control laws using two different configurations for controlling the nonlinear problem are investigated. As shown in §§ 4.1.2 and 4.1.3, the control kernels for full- and low-rank targets provide similar results, and we thus focus on the latter. We consider two scenarios, one using one sensor, one actuator and one target, and another using three of each. Although including statistics of the nonlinear terms within $\hat {{{\boldsymbol{\mathsf{F}}}}}$ in the nonlinear case would likely improve the results (Amaral et al. Reference Amaral, Cavalieri, Martini, Jordan and Towne2021), we chose not to pursue this so that the control law remains fixed as the forcing amplitudes increase, simplifying comparisons.

A series of nonlinear simulations were performed using spatially and temporally white forcing in the upstream region indicated in figure 7 with different forcing amplitudes defined by the root-mean-square (RMS) value $\epsilon$, i.e. $\hat {{{\boldsymbol{\mathsf{F}}}}} = \epsilon ^{2} {{\boldsymbol{\mathsf{I}}}}$. For reference, $\epsilon =10^{-2}$ leads to perturbations that reach 10–20 % of the baseflow velocity at the targets and, thus, is well into the nonlinear regime. As the control and estimation kernels are constructed to minimize the expected value of their respective cost functionals, a comparison should be based on an ensemble of simulations. Here we assume ergodicity of the system, and the ensemble averaging is replaced by time averaging. This assumption is validated using multiple runs for some configurations. The results for the linear problem presented in this section used an outflow condition on the rightmost domain, in order to properly compare the linear and nonlinear systems. The estimation and control laws should be based on a linear system that represents the nonlinear problem; however, for this problem, the impact of the different boundary conditions used in the linear and nonlinear systems is negligible, as reported next.

For small forcing amplitudes, the dynamics is dominated by linear mechanisms and, thus, it is expected that the performance of the estimation and control laws will be the same as demonstrated in the previous sections. Focusing initially on the control problem, representative snapshots of the flow and time evolution of perturbation norms for controlled and uncontrolled nonlinear systems are shown in figure 16. For the lower external forcing RMS, control of the nonlinear system is similar to the control of the linearized problem but, for larger amplitudes, it degrades. The perturbation norm, when normalized by the forcing RMS, reduces for larger forcing amplitudes due to nonlinear saturation. Estimation and control performances are measured as the non-estimated/controlled target energy fraction,

(4.2a,b)\begin{equation} \mathcal{E}_{est} = \frac{\displaystyle \sum_i \int (\boldsymbol{z}_i(t)-\boldsymbol{z}_{i,est}(t))^{2} \, {\rm d} t } {\displaystyle\sum_i \int (\boldsymbol{z}_i(t))^{2} \, {\rm d} t} ,\quad \mathcal{E}_{con} = \frac{\displaystyle \sum_i \int (\boldsymbol{z}_{i,con}(t))^{2} \, {\rm d} t}{\displaystyle\sum_i \int (\boldsymbol{z}_i(t))^{2} \, {\rm d} t}. \end{equation}

Figure 17 shows $\mathcal {E}_{est,con}$ as a function of $\epsilon$. Performances for the linear and the nonlinear systems are equivalent when small forcing amplitudes are used but, for the latter, it degrades for larger forcing amplitudes. The trend is also observed in figure 18, where time-series samples of a target for the uncontrolled problem, its causal estimation and for the controlled problem are shown. Comparing different runs, indicated by the black dots in figure 17, it can be observed that the spread is low for the low forcing amplitudes and increases for higher amplitudes, for which the control performance is degraded.

Figure 16.

(a–i) Snapshots of the nonlinear solutions for controlled and uncontrolled flows. Streamwise velocity is shown with the same colour scale as figure 7. Dotted/dashed contour lines indicate velocity excess/deficit of $\pm 5/\epsilon$ and $\pm 10/\epsilon$ from the baseflow. Circles, crosses and triangles indicate sensors, actuators and targets used. (j) Evolution of the perturbation norms, normalized by $\epsilon$.

Figure 17.

Control and estimation performances for the nonlinear (solid lines) and linearized (dashed) problems as a function of forcing amplitude. Black dots correspond to extra runs, with different random seeds for generating the external forcing. Results for the configuration shown in figure 16(i).

Figure 18.

Readings of the target located at $y=0$ for the controlled problem, it's estimate, and for the uncontrolled problem. The controlled problem corresponds to the configuration shown in figure 16(i). Results shown for (a) $\epsilon =10^{-4}$ and (b) $\epsilon =10^{-2}$.

Control and estimation performances for the linear problems are equivalent, indicating that all estimated perturbations are effectively controlled. A similar result is observed for the nonlinear problem disturbed by small external forcing. Estimation deteriorates for $\epsilon \approx 10^{-3}$, while control remains effective up to ${\epsilon \approx 10^{-2}}$. As the Kelvin–Helmholtz mode amplifies upstream disturbances, small but finite disturbances are amplified and can exhibit significant nonlinear dynamics. Such nonlinear effects, which are not accounted for by the linear approach used, degrade estimation performance. For the controlled system, these perturbations are cancelled before they are amplified, and, thus, the linear assumption is valid for a larger range of external forcing amplitudes for the controlled problem.

The degraded performance of the approaches presented here when applied to the nonlinear system disturbed with high-amplitude forcing is due to the saturation of the uncontrolled flow and to the violation of the assumption that the disturbances evolve linearly. From figure 16(j), comparing the normalized perturbation norm for the scenario with $\epsilon =10^{-2}$ to those with lower values of $\epsilon$, it can be seen that not only the perturbation norm for the controlled problem is larger, but also the norm of the uncontrolled problem is smaller. Both of these factors impact $\mathcal {E}_{con}$ as in (4.2a,b). The nonlinear interactions can be treated as additional external forcing terms in the linearized equation (McKeon & Sharma Reference McKeon and Sharma2010; Karban et al. Reference Karban, Bugeat, Martini, Towne, Cavalieri, Lesshafft, Agarwal, Jordan and Colonius2020). In this context, the degraded performance can be related to inadequacies of the forcing model. In particular, the nonlinear interactions are not, in general, zero mean when the linearization is performed around a baseflow, thus violating one of the assumptions used to construct the control laws. An alternative is to use a linearization around the mean flow, for which the nonlinear terms are, by definition, zero mean. In previous works (Martini et al. Reference Martini, Cavalieri, Jordan, Towne and Lesshafft2020; Amaral et al. Reference Amaral, Cavalieri, Martini, Jordan and Towne2021) it was shown that a more representative modelling of these forcing terms can lead to significant improvements in estimation, and, thus, potentially to improvements in control also for larger disturbances. Control formulations using a linearization around the mean flow may follow the methods in Högberg et al. (Reference Högberg, Bewley and Henningson2003) or Leclercq et al. (Reference Leclercq, Demourant, Poussot-Vassal and Sipp2019), with various mean flows used in successive linearisations of the system as control gains are progressively increased. Integrating these approaches into the control law is beyond the scope of this study, and will be considered in future work.