2.1. The resolvent operator

Let us consider a forced linear dynamical system

(2.1)\begin{equation} \frac{\mathrm{d}\kern0.7pt \boldsymbol{x}}{\mathrm{d} t}=\boldsymbol{A}\boldsymbol{x}+\boldsymbol{f}, \end{equation}

where $t\in \mathbb {R}$ denotes time, $\boldsymbol {x}(t)\in \mathbb {R}^n$ is the state whose dynamics are governed by the operator $\boldsymbol {A}\in \mathbb {R}^{n\times n}$, and $\boldsymbol {f}\in \mathbb {R}^n$ is the forcing. Such a system may arise from a semi-discretized partial differential equation, and in the case of fluid flows, the incompressible Navier–Stokes equations can be written in this form by projecting the velocity field onto a divergence-free basis to eliminate the pressure variable. The state $\boldsymbol {x}$ may represent either the deviation from a steady state of a laminar flow, or fluctuations about the temporal mean of a statistically stationary unsteady flow. In both cases, the matrix $\boldsymbol {A}$ is the linearization of the underlying nonlinear system about the corresponding base flow, either the equilibrium or mean flow. The forcing term may represent disturbances from the environment, model discrepancy, open-loop control actuation and/or the effect of nonlinear terms. Although we have defined $\boldsymbol {A}$, $\boldsymbol {x}$ and $\boldsymbol {f}$ to be real for simplicity, the methods presented below can be applied equally in the case of complex variables.

In the Laplace $s$-domain, the resolvent operator $\boldsymbol {H}(s)=(s\boldsymbol {I}-\boldsymbol {A})^{-1}$ is the transfer function from all possible inputs to all possible outputs in state space, with $\boldsymbol {I}\in \mathbb {R}^{n\times n}$ being the identity matrix. This operator encodes how any harmonic forcing $\hat {\boldsymbol {f}}(\omega )\,{\rm e}^{{\rm i}\omega t}$ at a specific frequency $\omega$ is amplified by the linear dynamics to produce a response

(2.2)\begin{equation} \hat{\boldsymbol{x}}(\omega)=\boldsymbol{H}({\rm i}\omega)\,\hat{\boldsymbol{f}}(\omega). \end{equation}

At any given frequency, a singular value decomposition (SVD) of the resolvent operator reveals the most responsive forcings, their gains, and the most receptive responses at that particular frequency. Specifically, the SVD factorizes the resolvent into $\boldsymbol {H}({\rm i}\omega )=\boldsymbol {\varPsi }(\omega )\,\boldsymbol {\varSigma }(\omega )\,\boldsymbol {\varPhi }(\omega )^*$, where $\boldsymbol {\varSigma } \in \mathbb {R}^{n\times n}$ is a diagonal matrix containing the resolvent gains $\sigma _1 \geqslant \sigma _2 \geqslant \cdots \geqslant \sigma _n \geqslant 0$, and $\boldsymbol {\varPhi }= [\boldsymbol {\phi }_1 \boldsymbol {\phi }_2\ \ldots \ \boldsymbol {\phi }_n]\in \mathbb {C}^{n\times n}$ and $\boldsymbol {\varPsi }=[\boldsymbol {\psi }_1 \boldsymbol {\psi }_2\ \ldots \ \boldsymbol {\psi }_n]\in \mathbb {C}^{n\times n}$ are unitary matrices whose columns are known as the forcing and response resolvent modes, $\boldsymbol {\phi }_j$ and $\boldsymbol {\psi }_j$, respectively.

2.2. The Gramians

Linear control theory deals with systems of the form

(2.3a,b)\begin{equation} \frac{\mathrm{d}\kern0.7pt \boldsymbol{x}}{\mathrm{d} t}=\boldsymbol{Ax}+\boldsymbol{Bu}, \quad \boldsymbol{y} = \boldsymbol{Cx},\end{equation}

where the matrices $\boldsymbol {B} \in \mathbb {R}^{n\times p}$ and $\boldsymbol {C}\in \mathbb {R}^{q \times n}$ are determined by the configuration of actuators and sensors, and describe how the control inputs $\boldsymbol {u}(t)\in \mathbb {R}^p$ act on the dynamics and which outputs $\boldsymbol {y}(t)\in \mathbb {R}^q$ are measured from the state, respectively. Note that $\boldsymbol {B}$, $\boldsymbol {C}$, $\boldsymbol {u}$ and $\boldsymbol {y}$ might be complex as well. In this setting, it is often important to identify the most controllable and observable directions in state space, that is, what are the states that are excited most efficiently by the inputs, and what are the states that are discerned most easily from the measured outputs. This is achieved by analysing the infinite-time observability and controllability Gramians, $\boldsymbol {W_o}$ and $\boldsymbol {W_c}$, defined as

(2.4a)$$\begin{gather} \boldsymbol{W_o} = \int_{0}^\infty {\rm e}^{\boldsymbol{A}^*t}\,\boldsymbol{C}^*\boldsymbol{C}\,{\rm e}^{\boldsymbol{A}t} \,\mathrm{d}t = \frac{1}{2{\rm \pi}}\int_{-\infty}^\infty \boldsymbol{H}({\rm i}\omega)^*\,\boldsymbol{C}^*\boldsymbol{C}\,\boldsymbol{H}({\rm i}\omega)\,\mathrm{d}\omega, \end{gather}$$

(2.4b)$$\begin{gather}\boldsymbol{W_c} = \int_{0}^\infty {\rm e}^{\boldsymbol{A}t}\,\boldsymbol{B}\boldsymbol{B}^*\,{\rm e}^{\boldsymbol{A}^*t} \,\mathrm{d}t = \frac{1}{2{\rm \pi}}\int_{-\infty}^\infty \boldsymbol{H}({\rm i}\omega)\,\boldsymbol{B}\boldsymbol{B}^*\,\boldsymbol{H}({\rm i}\omega)^*\, \mathrm{d}\omega, \end{gather}$$

where the rightmost expressions are the frequency domain representations which, as is well-known, feature the resolvent operator. It is important to point out that these representations are equivalent for stable systems only; whereas the integrals in the time domain diverge for the unstable case, the definitions in the frequency domain still hold as long as there are no marginal eigenvalues (Zhou et al. Reference Zhou, Salomon and Wu1999; Dergham et al. Reference Dergham, Sipp, Robinet and Barbagallo2011b). It is generally intractable to compute the Gramians directly using (2.4a) and (2.4b). Instead, they are often obtained by solving the Lyapunov equations

(2.5a)$$\begin{gather} \boldsymbol{A}^*\boldsymbol{W_o}+\boldsymbol{W_o A} + \boldsymbol{C}^*\boldsymbol{C}=\boldsymbol{0}, \end{gather}$$

(2.5b)$$\begin{gather}\boldsymbol{A W_c}+\boldsymbol{W_c} \boldsymbol{A}^* + \boldsymbol{B}\boldsymbol{B}^*=\boldsymbol{0}, \end{gather}$$

which can be achieved conveniently using readily available numerical routines (Skogestad & Postlethwaite Reference Skogestad and Postlethwaite2005; Brunton & Kutz Reference Brunton and Kutz2019). However, for high-dimensional systems, this becomes too computationally expensive, and typically the Gramians are approximated empirically from simulation data (Willcox & Peraire Reference Willcox and Peraire2002; Rowley Reference Rowley2005).

Both matrices, $\boldsymbol {W_o}$ and $\boldsymbol {W_c}$, are symmetric positive semi-definite, and their eigenvectors comprise hierarchically ordered orthogonal bases of known observable and controllable modes, respectively. For the case of full state outputs, the observability Gramian eigenvectors are the stochastic optimals (SOs) and, for the case of full state inputs, the eigenvectors of the controllability Gramian are the empirical orthogonal functions (EOFs) (Farrell & Ioannou Reference Farrell and Ioannou1993, Reference Farrell and Ioannou2001). Both SOs and EOFs also have deterministic interpretations that connect them to resolvent forcing and response modes. These connections are discussed in the remainder of this section.

2.3. SOs and EOFs are representative resolvent modes

We consider a linear system of the form in (2.1) by considering $\boldsymbol {C}=\boldsymbol {B}=\boldsymbol {I}$ such that we are measuring all possible outputs $\boldsymbol {y}=\boldsymbol {x}$ and forcing all possible inputs $\boldsymbol {u}=\boldsymbol {f}$. In this scenario, the Gramians become

(2.6a) $$\begin{gather} \boldsymbol{W_o} = \int_{0}^\infty {\rm e}^{\boldsymbol{A}{^*t}}\,{\rm e}^{\boldsymbol{A}t}\,\mathrm{d}t = \frac{1}{2{\rm \pi}}\int_{-\infty}^\infty \boldsymbol{H}({\rm i}\omega)^*\,\boldsymbol{H}({\rm i}\omega) \,\mathrm{d}\omega, \end{gather}$$

(2.6b)$$\begin{gather}\boldsymbol{W_c} = \int_{0}^\infty {\rm e}^{\boldsymbol{A}{t}}\,{\rm e}^{\boldsymbol{A}^*t}\,\mathrm{d}t = \frac{1}{2{\rm \pi}}\int_{-\infty}^\infty \boldsymbol{H}({\rm i}\omega)\,\boldsymbol{H}({\rm i}\omega)^* \,\mathrm{d}\omega. \end{gather}$$

These frequency domain expressions were leveraged by Dergham et al. (Reference Dergham, Sipp, Robinet and Barbagallo2011b) to obtain empirical approximations of the Gramians using data of long-time response to harmonic forcings. Furthermore, substituting the resolvent by its SVD $\boldsymbol {H}({\rm i}\omega )=\boldsymbol {\varPsi }(\omega )\,\boldsymbol {\varSigma }(\omega )\,\boldsymbol {\varPhi }(\omega )^*$, and using the orthogonality of the resolvent modes at each $\omega$, yields

(2.7a,b)\begin{equation} \boldsymbol{W_o} = \frac{1}{2{\rm \pi}}\int_{-\infty}^\infty \boldsymbol{\varPhi}(\omega)\,\boldsymbol{\varSigma}(\omega)^2\,\boldsymbol{\varPhi}(\omega)^* \, \mathrm{d}\omega, \quad \boldsymbol{W_c} = \frac{1}{2{\rm \pi}}\int_{-\infty}^\infty \boldsymbol{\varPsi}(\omega)\,\boldsymbol{\varSigma}(\omega)^2\,\boldsymbol{\varPsi}(\omega)^*\, \mathrm{d}\omega, \end{equation}

where no approximations have been made up to this point. We now consider a numerical quadrature to represent the integrals over the frequency domain as a sum over $m$ discrete frequencies $\omega _k$. For ease of notation, we consider a trapezoidal rule with a uniform frequency spacing $\Delta \omega$, but other quadratures and non-uniform spacing could be considered in general. Letting $\sigma _{jk}$, $\boldsymbol {\phi }_{jk}$ and $\boldsymbol {\psi }_{jk}$ denote the $j$th resolvent gain, forcing and response modes evaluated at the $k$th frequency leads to the following approximations of the Gramians:

(2.8a)\begin{align} \boldsymbol{W_o} &\approx \frac{\Delta \omega}{2{\rm \pi}}\sum_{k=1}^m \boldsymbol{\varPhi}(\omega_k)\,\boldsymbol{\varSigma}(\omega_k)^2\,\boldsymbol{\varPhi}(\omega_k)^* \end{align}

(2.8b)\begin{align} &=\frac{\Delta \omega}{2{\rm \pi}}\sum_{k=1}^m\left[\sigma_{1k}\boldsymbol{\phi}_{1k} \ \cdots \ \sigma_{nk}\boldsymbol{\phi}_{nk}\right] \left[\begin{array}{c} \sigma_{1k}\boldsymbol{\phi}_{1k}^* \\ \vdots\\ \sigma_{nk}\boldsymbol{\phi}_{nk}^* \end{array}\right] \end{align}

(2.8c)\begin{align} &= \frac{\Delta \omega}{2{\rm \pi}}\left[\sigma_{11}\boldsymbol{\phi}_{11} \ \cdots \ \sigma_{nm}\boldsymbol{\phi}_{nm}\right] \left[\begin{array}{c} \sigma_{11}\boldsymbol{\phi}_{11}^* \\ \vdots\\ \sigma_{nm}\boldsymbol{\phi}_{nm}^* \end{array}\right] \approx \boldsymbol{L_o}\boldsymbol{L_o}^* \end{align}

and

(2.9a)\begin{align} \boldsymbol{W_c} &\approx \frac{\Delta \omega}{2{\rm \pi}}\sum_{k=1}^m \boldsymbol{\varPsi}(\omega_k)\,\boldsymbol{\varSigma}(\omega_k)^2\,\boldsymbol{\varPsi}(\omega_k)^* \end{align}

(2.9b)\begin{align} &= \frac{\Delta \omega}{2{\rm \pi}}\sum_{k=1}^m\left[\sigma_{1k}\boldsymbol{\psi}_{1k} \ \cdots \ \sigma_{nk}\boldsymbol{\psi}_{nk}\right] \left[\begin{array}{c} \sigma_{1k}\boldsymbol{\psi}_{1k}^* \\ \vdots\\ \sigma_{nk}\boldsymbol{\psi}_{nk}^* \end{array}\right] \end{align}

(2.9c)\begin{align} &= \frac{\Delta \omega}{2{\rm \pi}}\left[\sigma_{11}\boldsymbol{\psi}_{11} \ \cdots \ \sigma_{nm}\boldsymbol{\psi}_{nm}\right] \left[\begin{array}{c} \sigma_{11}\boldsymbol{\psi}_{11}^* \\ \vdots\\ \sigma_{nm}\boldsymbol{\psi}_{nm}^* \end{array}\right] \approx \boldsymbol{L_c}\boldsymbol{L_c}^*, \end{align}

where the expressions in (2.8c) and (2.9c) reorder terms as the multiplication of two matrices, allowing us to recognize the Cholesky factors of the Gramians, $\boldsymbol {L_o}\approx (\Delta \omega /2{\rm \pi} )^{1/2}[\sigma _{11}\boldsymbol {\phi }_{11} \ \cdots \ \sigma _{nm}\boldsymbol {\phi }_{nm}]$ and $\boldsymbol {L_c}\approx (\Delta \omega /2{\rm \pi} )^{1/2}[\sigma _{11}\boldsymbol {\psi }_{11} \ \cdots \sigma _{nm}\boldsymbol {\psi }_{nm}]$. It is now clear that the column space of these Cholesky factors spans the same forcing or response subspace as all the corresponding resolvent modes at all frequencies. Moreover, because $\boldsymbol {L_o}$ and $\boldsymbol {L_c}$ are Cholesky factors, their left singular vectors are equivalent to the eigenvectors of $\boldsymbol {W_o}$ and $\boldsymbol {W_c}$. These eigenvectors are precisely the SOs and EOFs, and form orthogonal bases for these forcing and response subspaces. The framework introduced by Dergham et al. (Reference Dergham, Sipp and Robinet2013) uses the numerical quadrature described above, along with a rank truncation of the resolvent at every frequency, to approximate the leading SOs and EOFs for high-dimensional systems such as fluid flows.

The results presented above imply that SOs and EOFs can be interpreted as representative resolvent modes. Indeed, the leading SOs and EOFs will often closely resemble ad hoc selections of dominant – across frequencies – resolvent forcing and response modes, as shown in figure 1 for a turbulent channel flow example. Implementation details for this example can be found in § 5.

Figure 1. Forcing and response resolvent modes at handpicked high-gain frequencies compared to the leading SOs and EOFs for mean-flow-linearized minimal channel flow at $Re_\tau =185$. Real part of modes depicted as isosurfaces of wall-normal forcing and streamwise velocity response.

2.4. Gramian-based forcing and response modes

Although SOs and EOFs are usually explained as the structures that optimally contribute to and account for the sustained variance in stochastically excited flows, they also have an interesting deterministic interpretation (Farrell & Ioannou Reference Farrell and Ioannou2001). Specifically, the SOs and EOFs correspond to the most responsive forcings and the most receptive responses across all frequencies. This is formalized by deriving the SOs and EOFs as the solutions to the optimization problems presented below.

Let us first introduce the $\mathcal {L}_2$-norm for signals $\boldsymbol {y}(t)$. This norm measures energy integrated over time, as follows:

(2.10)\begin{equation} \|\boldsymbol{y}(t)\|_{\mathcal{L}_2}^2=\int_0^\infty \|\boldsymbol{y}(t)\|_2^2 \,\mathrm{d}t=\frac{1}{2{\rm \pi}}\int_{-\infty}^\infty \|\hat{\boldsymbol{y}}(\omega)\|_2^2\,\mathrm{d}\omega =\|\hat{\boldsymbol{y}}(\omega)\|_{\mathcal{L}_2}^2,\end{equation}

where $\hat {(\,\cdot\,)}$ denotes a Fourier-transformed variable, and the frequency domain expression follows from Plancherel's theorem.

In order to define our objective functions to be maximized, we consider the gain in the $\mathcal {L}_2$-norm for signals. Moreover, the last equality in (2.10) implies that, in the $\mathcal {L}_2$ sense, optimizing energy amplification integrated over all frequencies is equivalent to optimizing energy growth integrated over all time.

Hence the most responsive forcing across all frequencies is the mode $\boldsymbol {v}$ that, mapped through the resolvent, produces the most amplified response in the $\mathcal {L}_2$ sense and is found by maximizing the gain

(2.11)\begin{equation} \frac{\|\boldsymbol{H}({\rm i}\omega)\,\boldsymbol{v} \|^2_{\mathcal{L}_2}}{\|\boldsymbol{v}\|_2^2} = \frac{\|{\rm e}^{\boldsymbol{A}t}\,\boldsymbol{v} \|^2_{\mathcal{L}_2}}{\|\boldsymbol{v}\|_2^2} = \frac{\boldsymbol{v}^* \displaystyle\int\nolimits_0^\infty {\rm e}^{\boldsymbol{A}^*t}\,{\rm e}^{\boldsymbol{A}t} \,\mathrm{d}t \,\boldsymbol{v}}{\boldsymbol{v}^*\boldsymbol{v}} = \frac{\boldsymbol{v}^*\boldsymbol{W_o} \boldsymbol{v}}{\boldsymbol{v}^*\boldsymbol{v}}, \end{equation}

where the observability Gramian $\boldsymbol {W_o}$ emerges naturally from expanding the norm definition. In fact, the cost function in the rightmost expression corresponds to a Rayleigh quotient, thus the optimizer is the leading eigenvector of $\boldsymbol {W_o}$, that is, the leading SO. Subsequent eigenvectors then solve for the next most responsive forcings. Although the association between the most observable modes and the most responsive forcings may seem counter-intuitive, these are the structures that generate the most energetic responses, and can therefore be most easily discerned from the time history of the state.

Now, the most receptive response across all frequencies is the mode $\boldsymbol {v}$ that best aligns, in the $\mathcal {L}_2$ sense, with all possible responses, which are spanned by the columns of the resolvent, thus maximizing the gain

(2.12)\begin{equation} \frac{\|\boldsymbol{H}({\rm i}\omega)^*\,\boldsymbol{v} \|^2_{\mathcal{L}_2}}{\|\boldsymbol{v}\|_2^2} = \frac{\|{\rm e}^{\boldsymbol{A}^*t}\,\boldsymbol{v} \|^2_{\mathcal{L}_2}}{\|\boldsymbol{v}\|_2^2} = \frac{\boldsymbol{v}^* \displaystyle\int\nolimits_0^\infty {\rm e}^{\boldsymbol{A}t}\,{\rm e}^{\boldsymbol{A}^*t} \,\mathrm{d}t \,\boldsymbol{v}}{\boldsymbol{v}^*\boldsymbol{v}} = \frac{\boldsymbol{v}^*\boldsymbol{W_c} \boldsymbol{v}}{\boldsymbol{v}^*\boldsymbol{v}}, \end{equation}

where now a Rayleigh quotient with the controllability Gramian $\boldsymbol {W_c}$ appears in the last expression. Thus the optimizer corresponds to its leading eigenvector, that is, the leading EOF. Again, subsequent eigenvectors then solve for the next most receptive responses in the same sense. The intuition behind the association between controllable modes and the most receptive responses lies in the fact that these are the states that the system is most easily steered to through the action of any input or forcing.

When working with modal decompositions of fluid flows (or dynamical systems in general), it is useful to think about the optimization problems being solved to provide physical interpretations for the modes. This has been particularly important for non-normal systems, where non-modal stability analyses have provided valuable insights regarding physical mechanisms for linear energy amplification (Schmid Reference Schmid2007). As we have shown, and summarize in table 1, the forcing and response modes produced by the Gramian eigenvectors are closely connected to those obtained from transient growth and resolvent analyses.

Table 1. Forcing and response modes for non-normal systems, and the optimization problems that they solve.

It is important to note that we have been using the Euclidean $2$-norm to represent the energy of a state vector $\|\boldsymbol {v}\|_2^2$. It is often the case that a physically meaningful inner product considers a positive-definite weighting matrix $\boldsymbol {Q}$ so that the relevant metric is $\|\boldsymbol {v}\|_{\boldsymbol {Q}}^2=\boldsymbol {v}^*\boldsymbol {Qv}$. This weighting might account for integration quadratures or scaling of heterogeneous variables in multi-physics systems (Herrmann-Priesnitz et al. Reference Herrmann-Priesnitz, Calderón-Muñoz, Diaz and Soto2018a). However, in this scenario, one can modify the problem conveniently and easily to work with the Euclidean $2$-norm by taking the state as $\boldsymbol {Fv}$ and the dynamics as $\boldsymbol {F A F}^{-1}$, where $\boldsymbol {F}^*$ is the Cholesky factor of $\boldsymbol {Q}=\boldsymbol {F}^*\boldsymbol {F}$.

3.1. Orthogonal input and output projections

The goal of input and output projections is to find projectors $\mathbb {P}$ to produce low-rank approximations of either inputs or outputs, respectively, such that the input–output behaviour remains as close as possible to that of the original system of form (2.3a,b). In the case of input projections, the approximated dynamics is governed by

(3.1a,b)\begin{equation} \frac{\mathrm{d} \tilde{\boldsymbol{x}}}{\mathrm{d} t}=\boldsymbol{A}\tilde{\boldsymbol{x}}+\boldsymbol{B}\mathbb{P}\boldsymbol{u}, \quad \tilde{\boldsymbol{y}} = \boldsymbol{C}\tilde{\boldsymbol{x}}, \end{equation}

and for the case of output projections, it is described by

(3.2a,b)\begin{equation} \frac{\mathrm{d}\kern0.7pt \boldsymbol{x}}{\mathrm{d} t}=\boldsymbol{Ax}+\boldsymbol{Bu}, \quad \tilde{\boldsymbol{y}} = \mathbb{P}\boldsymbol{Cx}. \end{equation}

In both scenarios, we want the outputs $\boldsymbol {y}(t)$ and $\tilde {\boldsymbol {y}}(t)$ for the original and projected systems to be close for a sequence of inputs $\boldsymbol {u}(t)$. We will consider the $\mathcal {L}_2$-norm for signals, introduced in (2.10), as our metric of choice, which measures energy integrated over time. Under this norm, the input projection error becomes

(3.3a)\begin{align} \|\boldsymbol{y}(t)-\tilde{\boldsymbol{y}}(t)\|_{\mathcal{L}_2}^2&= \frac{1}{2{\rm \pi}}\int_{-\infty}^\infty \|\boldsymbol{G}({\rm i}\omega)\hat{\boldsymbol{u}}(\omega)-\boldsymbol{G}({\rm i}\omega)\,\mathbb{P}\,\hat{\boldsymbol{u}}(\omega)\|_2^2\, \mathrm{d}\omega \end{align}

(3.3b)\begin{align} &= \frac{1}{2{\rm \pi}}\int_{-\infty}^\infty \|\boldsymbol{G}({\rm i}\omega)\,(\boldsymbol{I}-\mathbb{P})\,\hat{\boldsymbol{u}}(\omega)\|_2^2\,\mathrm{d}\omega, \end{align}

and the output projection error becomes

(3.4a)\begin{align} \|\boldsymbol{y}(t)-\tilde{\boldsymbol{y}}(t)\|_{\mathcal{L}_2}^2 &= \frac{1}{2{\rm \pi}}\int_{-\infty}^\infty \|\boldsymbol{G}({\rm i}\omega)\,\hat{\boldsymbol{u}}(\omega)-\mathbb{P}\,\boldsymbol{G}({\rm i}\omega)\,\hat{\boldsymbol{f}}(\omega)\|_2^2\,\mathrm{d}\omega \end{align}

(3.4b)\begin{align} &= \frac{1}{2{\rm \pi}}\int_{-\infty}^\infty \|(\boldsymbol{I}-\mathbb{P})\,\boldsymbol{G}({\rm i}\omega)\,\hat{\boldsymbol{u}}(\omega)\|_2^2\,\mathrm{d}\omega, \end{align}

where $\boldsymbol {G}(s)=\boldsymbol {C}(s\boldsymbol {I}-\boldsymbol {A})^{-1}\boldsymbol {B}=\boldsymbol {C}\,\boldsymbol {H}(s)\,\boldsymbol {B}$ is the transfer function for the unprojected system (2.3a,b).

The values of $\mathbb {P}$ that minimize the errors (3.3) and (3.4) depend on the class of inputs admitted and on the type of projector sought. The methods introduced by Dergham et al. (Reference Dergham, Sipp and Robinet2011a) and Rowley (Reference Rowley2005) seek orthogonal projectors to minimize the errors (3.3) and (3.4), respectively, for inputs that are impulses on every possible degree of freedom. Given an orthonormal set of vectors $\boldsymbol {V}=[\boldsymbol {v}_{1} \ \cdots \ \boldsymbol {v}_{r}]$, an orthogonal projector onto the column space of $\boldsymbol {V}$ is given by $\mathbb {P}_{\perp }=\boldsymbol {VV}^*$.

If such a rank-$r$ projector minimizes the input projection error (3.3), then the modes in $\boldsymbol {V}$ correspond to the most dynamically relevant inputs, that is, the response of the system can be well predicted when approximating $\boldsymbol {u}(t)\approx \mathbb {P}_{\perp }\,\boldsymbol {u}(t)= \boldsymbol {VV}^*\,\boldsymbol {u}(t)=\boldsymbol {V}\,\boldsymbol {d}(t)$, with $\boldsymbol {d}(t)\in \mathbb {C}^{r}$. In addition to the physical insight regarding the most responsive inputs, these modes can be leveraged to design spatially distributed actuators. Building actuators that are tailored to excite these modes results in a controller that can act along the most responsive directions in state space. The optimal basis for orthogonal input projection is given by the leading POD modes from adjoint impulse response simulations, which, in the case of full state inputs, are equivalent to the leading observable modes (Dergham et al. Reference Dergham, Sipp and Robinet2011a).

Similarly, if the output projection error in (3.4) is minimized, then $\boldsymbol {V}$ spans the subspace of the most receptive measurements. In this case, the modes in $\boldsymbol {V}$ allow for a low-dimensional representation of the output in terms of mode amplitudes, that is, $\boldsymbol {y}(t)\approx \mathbb {P}_{\perp }\,\boldsymbol {y}(t)= \boldsymbol {VV}^*\,\boldsymbol {y}(t)=\boldsymbol {V}\,\boldsymbol {a}(t)$, with $\boldsymbol {a}(t)\in \mathbb {C}^{r}$. The optimal basis for orthogonal output projection is given by the leading POD modes from direct impulse response simulations, which, in the case of full state outputs, are equivalent to the leading controllable modes (Rowley Reference Rowley2005; Holmes et al. Reference Holmes, Lumley, Berkooz and Rowley2012).

Furthermore, in this work we are concerned with fluid flows that can be forced anywhere and measured everywhere, that is, the system has full state inputs and outputs simultaneously. Recall that if this is the case, then the transfer function becomes the resolvent operator $\boldsymbol {G}(s)=\boldsymbol {H}(s)$. In this scenario, the optimal bases that minimize (3.3) and (3.4) are comprised of flow structures that coincide with the leading SOs and EOFs, and can be interpreted as the most responsive forcings and receptive responses, respectively.

3.2. Interpolatory projectors

We have presented the idea of using projectors to produce low-rank approximations of inputs and outputs, or forcings and responses in dynamical systems of the form in (2.1) with full state inputs and outputs. The optimal orthogonal projectors that minimize the errors (3.3) and (3.4) in that case use the leading SOs and EOFs as bases. Here, we define another kind of projector – an interpolatory projector – that, as we show later, is useful practically for sensor and actuator placement and open-loop control.

An interpolatory projector is an oblique projector with the special property of preserving exactly certain entries in any vector upon which it acts. These entries are referred to as interpolation points, thus justifying the name, which was first introduced in Chaturantabut & Sorensen (Reference Chaturantabut and Sorensen2010) and later adopted by more recent work (Sorensen & Embree Reference Sorensen and Embree2016; Gidisu & Hochstenbach Reference Gidisu and Hochstenbach2021).

Definition 3.1 Given a full-rank matrix $\boldsymbol {V}\in \mathbb {R}^{n\times r}$ and a set of distinct indices $\boldsymbol {\gamma }\in \mathbb {N}^r$, the interpolatory projector for $\boldsymbol {\gamma }$ onto the column space of $\boldsymbol {V}$ is

(3.5)\begin{equation} \mathbb{P}_{:}\equiv \boldsymbol{V}(\boldsymbol{P}^{\rm{T}} \boldsymbol{V})^{{-}1}\boldsymbol{P}^{\rm{T}}, \end{equation}

provided that $\boldsymbol {P}^{\rm {T}} \boldsymbol {V}$ is invertible, where $\boldsymbol {P}=[\boldsymbol {e}_{{\gamma }_1} \ \cdots \boldsymbol {e}_{{\gamma }_r}] \in \mathbb {R}^{n\times r}$ contains $r$ columns of the identity matrix given by $\boldsymbol {\gamma }$.

The matrix $\boldsymbol {P}$ is a sparse matrix, referred to as the sampling matrix because it samples the entries given by the indices $\boldsymbol {\gamma }$ of any vector $\boldsymbol {x}\in \mathbb {R}^n$ upon which it acts, that is, $\boldsymbol {P}^{\rm {T}}\boldsymbol {x}=\boldsymbol {x}(\boldsymbol {\gamma })$. Just like an orthogonal projector $\mathbb {P}_{\perp }=\boldsymbol {VV}^{\rm {T}}$, when $\mathbb {P}_{:}$ is applied to a vector, the result is in the column space of $\boldsymbol {V}$. However, $\mathbb {P}_{:}$ is an oblique projector, and it also has the following distinct property: for any $\boldsymbol {x}\in \mathbb {R}^n$,

(3.6)\begin{equation} (\mathbb{P}_{:}\boldsymbol{x})(\boldsymbol{\gamma})=\boldsymbol{P}^{\rm{T}}\mathbb{P}_{:}\boldsymbol{x}=\boldsymbol{P}^{\rm{T}}\boldsymbol{V}(\boldsymbol{P}^{\rm{T}} \boldsymbol{V})^{{-}1} \boldsymbol{P}^{\rm{T}}\boldsymbol{x} = \boldsymbol{P}^{\rm{T}}\boldsymbol{x} = \boldsymbol{x}(\boldsymbol{\gamma}), \end{equation}

so the projected vector $\mathbb {P}_{:}\boldsymbol {x}$ matches $\boldsymbol {x}$ in the entries indicated by $\boldsymbol {\gamma }$. These entries are known as the interpolation or sampling points, and are determined by the non-zero entries in the sampling matrix. A simple example of interpolatory projection is presented in Appendix A.

3.3. Proposed method

For our interpolatory input and output projections, we take the leading SOs and EOFs as the basis $V$, respectively. It is important to note that these projectors are sub-optimal with respect to (3.3) and (3.4), and therefore will never show better performance (under these metrics) than the orthogonal projectors of the same rank built with the same bases. Furthermore, the choice of the sampling locations, or equivalently the design of $\boldsymbol {P}$, is critical for the performance of the interpolatory projection to approximate that of the orthogonal one. In addition, when using the leading SOs or EOFs in $\boldsymbol {V}$, the interpolation points that result in small projection errors are physically interesting, since they represent responsive or receptive spatial locations, respectively. However, finding the optimal sampling matrix that minimizes the projection error for a given basis $\boldsymbol {V}$ requires a brute force search over all possible sampling location combinations. While this can be achieved for small-scale systems (Chen & Rowley Reference Chen and Rowley2011), with $n\sim {O}(10^4)$ at the time of writing, the computational cost makes it intractable in higher dimensions. Alternatively, there are several fast greedy algorithms that can be used to approximate the optimal sensor locations, avoiding the combinatorial search. In the context of reduced-order modelling, the empirical and discrete empirical interpolation methods, EIM (Barrault et al. Reference Barrault, Maday, Nguyen and Patera2004) and DEIM (Chaturantabut & Sorensen Reference Chaturantabut and Sorensen2010), were developed to find locations to interpolate nonlinear terms in a high-dimensional dynamical system, which is known as hyper-reduction. An even simpler and equally efficient approach is the Q-DEIM algorithm, introduced by Drmac & Gugercin (Reference Drmac and Gugercin2016), that leverages the pivoted QR factorization to select the sampling points. Manohar et al. (Reference Manohar, Brunton, Kutz and Brunton2018) showed that this is also a robust strategy for sparse sensor placement for state reconstruction on a range of applications. Here, we adhere to this method because it provides near-optimal interpolation points and is simple to implement by leveraging the pivoted QR factorization available in most scientific computing software packages. The pseudo-code for this greedy sensor selection strategy, adapted from Manohar et al. (Reference Manohar, Brunton, Kutz and Brunton2018), is shown in Algorithm 1, where we use Matlab's array slicing notation.

Algorithm 1 Greedy sensor selection using pivoted QR (Manohar et al. 2018)

Therefore, the approach used to build $\boldsymbol {P}$ for our interpolatory input or output projections consists in taking a pivoted QR factorization of the transpose of the basis $V$, which is comprised of the leading either SOs or EOFs, respectively. Then the first $r$ QR pivots correspond to the sampling points, or sensor locations, which we will refer to as tailored sensors since they depend solely on the corresponding $\boldsymbol {V}$, and are therefore tailored specifically to that basis (Manohar et al. Reference Manohar, Brunton, Kutz and Brunton2018).

Combining the SOs and EOFs with the pivoted QR factorization for sampling point selection, we obtain an equation-driven framework for interpolatory input and output projections that is easy to implement. In the following, we explain how this approach can be leveraged for several control-related tasks.

3.4. Sparse sensor and actuator placement

The leading SOs span the subspace of the dominant forcings for response prediction. Using an interpolatory projector, these dominant forcings can be approximated as $\boldsymbol {f}(t)\approx \mathbb {P}_{:}\,\boldsymbol {f}(t)=\boldsymbol {V}(\boldsymbol {P}^{\rm {T}} \boldsymbol {V})^{-1} \boldsymbol {P}^{\rm {T}}\,\boldsymbol {f}(t)=\boldsymbol {V}(\boldsymbol {P}^{\rm {T}} \boldsymbol {V})^{-1}\,\boldsymbol {d}(t)$, where now $\boldsymbol {d}$ corresponds to the entries of $\boldsymbol {f}$ at the $r$ interpolation points. This means that we can approximate the dynamically relevant parts of the forcing from sparse sensor measurements, and if a model of the dynamics is available, predict the response of the system. Moreover, the interpolation points also make attractive open-loop actuator locations to act on the SOs and mitigate the disturbances with the largest impact on the dynamics. Similarly, using the EOFs, we can reconstruct the response $\boldsymbol {x}(t)\approx \mathbb {P}_{:}\,\boldsymbol {x}(t)=\boldsymbol {V}(\boldsymbol {P}^{\rm {T}} \boldsymbol {V})^{-1} \boldsymbol {P}^{\rm {T}}\,\boldsymbol {x}(t)=\boldsymbol {V}(\boldsymbol {P}^{\rm {T}} \boldsymbol {V})^{-1}\,\boldsymbol {a}(t)$ from sparse sensor measurements $\boldsymbol {a}$ corresponding to the entries of $\boldsymbol {x}$ at the $r$ interpolation points.

Balanced modes, corresponding to the most jointly controllable and observable states of a linear dynamical system (Moore Reference Moore1981; Rowley & Dawson Reference Rowley and Dawson2017), have also been leveraged for sensor and actuator placement by Manohar et al. (Reference Manohar, Kutz and Brunton2021). However, that work was focused on feedback control, whereas here we are proposing sensors for forcing or response reconstruction and actuators for open loop control. Although balanced modes may also seem like attractive choices of bases to build input and output projectors, the subspaces that they span are sub-optimal for projection in terms of the errors (3.3) and (3.4), since the balancing implies a trade-off between observability and controllability. Therefore, for the tasks presented in this work, that rely on input and output projections, the eigenvectors of the Gramians are more principled choices of bases. This way, the resulting interpolatory projectors are approximating the optimal orthogonal projectors.

3.5. Partial disturbance feedforward control

Disturbance feedforward control is an open-loop control strategy that uses measurements of the disturbances to a system and a model of its dynamics to determine a sequence of control inputs with the objective of mitigating the effect of disturbances. This strategy is used typically to reject low-frequency, high-amplitude disturbances and perform aggressive reference tracking, in combination with a feedback controller that attenuates high-frequency sensor noise and compensates for model uncertainty (Bade et al. Reference Bade, Hanson, Belson, Naguib, Lavoie and Rowley2016; Papadakis, Lu & Ricco Reference Papadakis, Lu and Ricco2016; Williams & King Reference Williams and King2018; Tol, Kotsonis & De Visser Reference Tol, Kotsonis and De Visser2019; Herrmann et al. Reference Herrmann, Brunton, Pohl and Semaan2022). Consider a linear dynamical system

(3.7a,b)\begin{equation} \frac{\mathrm{d}\kern0.7pt \boldsymbol{x}}{\mathrm{d} t}=\boldsymbol{Ax}+\boldsymbol{w}+\boldsymbol{Bu}, \quad \boldsymbol{y}=\boldsymbol{Cw},\end{equation}

that is forced simultaneously by exogenous disturbances $\boldsymbol {w}$ and the control actuation term $\boldsymbol {Bu}$, with $\boldsymbol {u}\in \mathbb {C}^p$ and $\boldsymbol {B}\in \mathbb {C}^{n\times p}$, where $\boldsymbol {y}=\boldsymbol {Cw}$ are measurements of the disturbances through the linear mapping $\boldsymbol {C}\in \mathbb {C}^{q\times n}$. This is a particular case of the disturbance feedforward control problem presented in Doyle et al. (Reference Doyle, Glover, Khargonekar and Francis1988) with no feedback, where the error signal $\boldsymbol {z}=\boldsymbol {x}$ is just the state. Therefore, for a given set-up, an optimal control law can be derived within the framework of $\mathcal {H}_2$ control (Doyle et al. Reference Doyle, Glover, Khargonekar and Francis1988; Skogestad & Postlethwaite Reference Skogestad and Postlethwaite2005). However, we consider the joint problem of designing simultaneously the actuation matrix $\boldsymbol {B}$, the measurement matrix $\boldsymbol {C}$ and a simple control law with a constant gain matrix $\boldsymbol {u}=-\boldsymbol {Ky}$, with the aim of mitigating only the most meaningful disturbances. Substituting the control law in (3.7a,b), the controlled system dynamics becomes

(3.8a,b)\begin{equation} \frac{\mathrm{d}\kern0.7pt \boldsymbol{x}}{\mathrm{d} t}=\boldsymbol{Ax}+\boldsymbol{f}, \quad \boldsymbol{f} = \left(\boldsymbol{I}-\boldsymbol{BKC}\right)\boldsymbol{w},\end{equation}

where the system is forced by the effective forcing $\boldsymbol {f}$. Now we want to find a combination of actuators $\boldsymbol {B}$, disturbance sensors $\boldsymbol {C}$ and gain $\boldsymbol {K}$ that minimizes the effect of $\boldsymbol {f}$ on the response of the system. This is achieved if $(\boldsymbol {I}-\boldsymbol {BKC})$ projects $\boldsymbol {w}$ onto the least responsive disturbances, or equivalently, if $\boldsymbol {BKC}$ projects onto the most responsive disturbances, given by the leading SOs of the system (3.8a,b). In other words, the sensors, control gains and actuators can be selected conveniently so that the matrix product of $\boldsymbol {B}$, $\boldsymbol {K}$ and $\boldsymbol {C}$ produces an orthogonal or interpolatory input projector. Note that this approach only partially mitigates the disturbances, which is why we refer to it as partial disturbance feedforward control.

This approach is particularly attractive for noise-amplifier systems, where a few preferred forcing structures are responsible for most of the energy sustained by the response (Dergham et al. Reference Dergham, Sipp and Robinet2011a, Reference Dergham, Sipp and Robinet2013). In § 4, we demonstrate the performance of partial disturbance feedforward control on a stochastically forced linear system considering different types of sensors and actuators.

Furthermore, the proposed approach may also prove to be interesting in the context of turbulent flows. Specifically, for wall-bounded turbulent flows, the SOs of the mean-flow linearization are expected to be an efficient basis for input projection because of their behaviour as selective noise amplifiers (Dergham et al. Reference Dergham, Sipp and Robinet2013). However, the effect of control inputs would cause mean-flow deformation. Whether or not the actuators designed based on the SOs of the uncontrolled flow would remain effective when control is in place has yet to be investigated. This will be the subject of future research and is outside the scope of the present work.

4.1. Most responsive disturbances

For the set of parameters used, the linearized complex Ginzburg–Landau equation is known to behave as a selective amplifier of disturbances (Bagheri et al. Reference Bagheri, Henningson, Hoepffner and Schmid2009b). When forced with white noise, the dynamics of the system filters disturbances and generates a response where most of the energy is contained in coherent structures, as shown in figure 2(a). This implies that, among all possible ways to excite the system, there are a few preferred forcing patterns that are dynamically relevant and have a dominant effect on the produced response.

Figure 2. (a) Response of the linearized complex Ginzburg–Landau equation to white noise forcing. The system behaves as a selective disturbance amplifier. (b) Two bases of representative forcing modes are used to build orthogonal input projectors: the leading forcing POD modes computed from snapshots of the disturbances, and the leading SOs. (c) Projected forcings and their induced responses using $20$ forcing POD modes versus $3$ SOs. (d) Normalized input projection error as a function of the projector rank $r$. The real part is displayed for forcing, response and mode plots. These findings follow directly from the work of Dergham et al. (Reference Dergham, Sipp and Robinet2011a).

We build orthogonal projectors onto two different subspaces, spanned by the bases $\boldsymbol {V}$ shown in figure 2(b), that one might consider when seeking low-rank approximations of the forcing. For the first subspace, we gather snapshots of the forcing, perform POD on those snapshots, and retain the leading $r$ forcing POD modes as our basis to build a rank-$r$ projector. The second subspace is the one spanned by the leading $r$ eigenvectors of the observability Gramian. The projection of the forcing onto these subspaces, and the response of the system to the projected forcings, are shown in figure 2(c) for $r=20$ POD modes and $r=3$ $\boldsymbol {W_o}$ eigenmodes.

The projected forcing using POD resembles closely the original forcing, which is expected since POD provides the optimal low-rank approximation based on empirical observations. Despite this, the response is predicted more accurately using the input projection onto the SOs, even when using just $3$ modes compared to the $20$ used for the POD-based projection, as shown in figure 2(c). More quantitatively, the error in the response prediction, measured by the normalized $\mathcal {L}_2$-norm, is always lower for the same forcing projection rank $r$, and decreases much more rapidly for low-rank approximations, as shown in figure 2(d).

This performance gap is expected since the POD modes of the disturbance field, albeit optimal at reconstructing the observed disturbances, are completely agnostic to the dynamics of the system. The purpose of the comparison is to highlight the dynamical significance of the SOs as the most responsive disturbances, providing a subspace for input projection that is optimal for response prediction. We remark that the findings presented in this subsection are not new since they follow directly from the application of the method developed by Dergham et al. (Reference Dergham, Sipp and Robinet2011a).

4.2. Sensor placement for prediction

Using interpolatory input projection, the most responsive disturbances buried within white noise forcing can be exposed successfully from sparse sensor measurements, as shown in figure 3(a) for the case of $8$ sensors. The projected forcing captures the dynamically relevant parts of the disturbances, allowing for an accurate prediction of the energetic coherent structures in the response, as also shown in figure 3(a). Moreover, the error induced by the interpolatory projection on the response prediction rapidly approaches the error induced by an orthogonal projector, built using the same basis, as the basis size increases, as shown in figure 3(b). The spatial locations found for the sensors, tailored to the eigenvectors of $\boldsymbol {W_o}$, resemble closely the optimal placements computed by Chen & Rowley (Reference Chen and Rowley2011). However, our approach circumvents solving the computationally expensive optimization problem presented in that work by using the pivoted QR method. Furthermore, the sensor locations found are also similar to those found by Manohar et al. (Reference Manohar, Kutz and Brunton2021), which is expected since those were tailored to the adjoint balanced modes, which span almost the same subspace as the eigenvectors of the observability Gramian for this particular case. Importantly, in that and other previous works, the sensor placement was optimized for feedback control, whereas here, sensors aim to reconstruct the disturbance field to use it for prediction. As discussed in § 3, for this purpose, balanced modes are sub-optimal and SOs are the optimal choice.

Figure 3. (a) Forcing applied to the linearized complex Ginzburg–Landau equation and its induced response compared to an interpolatory projection of the forcing and its induced response. The projector interpolates the forcing in the span of the leading $r=8$ SOs from measurements at $r=8$ spatial locations (sensors). (b) Normalized input projection error using interpolatory (dashed) and orthogonal (solid) projectors as a function of the projector rank $r$. (c) Sensor locations as a function of the size of the basis used for projection. Sampling points are selected using pivoted QR, as in Manohar et al. (Reference Manohar, Brunton, Kutz and Brunton2018), and are therefore tailored to the basis of SOs. Plots show the real parts of forcings and responses.

4.3. Open-loop control design

Here, we consider the application of the partial disturbance feedforward control approach described in § 3 to the linearized complex Ginzburg–Landau equation subject to white noise disturbances. This strategy leverages the concept input projections to address only the most meaningful disturbances. With this strategy in mind, we design three feedforward controllers considering different restrictions on the allowable sensors and actuators.

First, we assume that we have access to full disturbance measurements, that is, $\boldsymbol {C}=\boldsymbol {I}$, and that spatially distributed (body forcing) actuators are feasible. We design our actuators to excite directly the most responsive disturbances, that is, $\boldsymbol {B}=\boldsymbol {V}$, with $\boldsymbol {V}$ containing the leading SOs. Therefore, choosing the gain to be $\boldsymbol {K}=\boldsymbol {V}^*$ results in $\boldsymbol {BKC}$ being the optimal orthogonal input projector. Considering $r=10$ distributed actuators, the feedforward controller is able to successfully reject disturbances by projecting-out the components that would cause the most energetic response, as shown in figure 4(b). More precisely, a $97\,\%$ reduction in the $\mathcal {L}_2$-norm of the fluctuations is achieved.

Figure 4. (a) Feedforward control strategy to reject disturbances in the linearized complex Ginzburg–Landau equation. The control actuation and the response of the controlled system are shown for three sensor and actuator configurations: (b) full disturbance measurements and $10$ spatially distributed (body forcing) actuators; (c) $10$ point sensors and $10$ spatially distributed actuators; and (d) $10$ localized sensors and actuators. Plots for the forcing, control actuations and responses display the real parts of these fields.

For a second case, we allow only pointwise measurements of the disturbances, that is, $\boldsymbol {C}=\boldsymbol {P}^{\rm {T}}$, and still consider spatially distributed actuators. Choosing our actuators to be $\boldsymbol {B}=\boldsymbol {V}$ means that the gain $\boldsymbol {K}=(\boldsymbol {P}^{\rm {T}}\boldsymbol {V})^{-1}$ makes $\boldsymbol {BKC}$ an interpolatory input projector. In this scenario, the full disturbance field is interpolated from sparse sensors, and then the actuators attempt to cancel out the components that would lead to the most energetic response. Using $r=10$ spatially distributed actuators and $10$ point sensors, selected via pivoted QR, the feedforward controller has only slightly lower performance than when using full measurements, as shown in figure 4(c). Specifically, a $93\,\%$ reduction in the $\mathcal {L}_2$-norm of the fluctuations is obtained in this case.

Finally, we examine a more realistic scenario where only spatially localized sensors and actuators are allowed. Specifically, we consider $\boldsymbol {B}=\tilde {\boldsymbol {P}}$ and $\boldsymbol {C}=\tilde {\boldsymbol {P}}^{\rm {T}}$, where, rather than $r$ columns of the identity as in $\boldsymbol {P}$, the columns of $\tilde {\boldsymbol {P}}$ are given by discretized Gaussian functions

(4.2)\begin{equation} \exp\left( -\frac{\left(x-x_j \right)^2}{2\sigma^2}\right), \quad j=1, \dots, r, \end{equation}

localized at sampling points $x_j$, where we choose $\sigma =0.6$. Using a control gain $\boldsymbol {K}=(\tilde {\boldsymbol {P}}^{\rm {T}}\tilde {\boldsymbol {P}})^{-1}$, the actuation then attempts to cancel out directly the disturbances measured in the regions localized around the sampling points. We hypothesize that the optimal sensor locations for response prediction are also good for actuator placement, since they correspond to the optimal locations to interpolate the most responsive forcings. Therefore, leveraging the first $r=10$ SOs, we use pivoted QR to select $10$ sampling points to place disturbance sensors and actuators. The resulting controller has a reasonable disturbance rejection performance, as shown in figure 4(d). In this case, an $89\,\%$ reduction in the $\mathcal {L}_2$-norm of the fluctuations is achieved.

5.1. Most receptive perturbations

We examine the suitability of the EOFs and POD modes as bases $\boldsymbol {V}$ for the dominant response subspace. We build rank-$r$ orthogonal projectors from the leading $r$ elements in these bases, and use them to project the DNS snapshots. Qualitatively, both POD modes and $\boldsymbol {W_c}$ eigenmodes seem to capture the main coarse features in the flow with a rank $100$ projection and further details as $r$ increases, as shown in figure 5(a) for the streamwise velocity.

Second-order statistics of the velocity fluctuations are approached by those obtained for the projected snapshots as the rank of the projector increases with both bases, as shown in figure 5(c). However, there is a slower convergence for the Gramian-based projection. Specifically, although the Gramian eigenmodes capture reasonably the peak in the root-mean-square (r.m.s.) profiles, a higher number of modes is required to resolve properly the statistics close to the wall and near the channel centre.

More quantitatively, the response projection error, measured by the normalized $\mathcal {L}_2$-norm, drops faster with increasing $r$ and is always lower when using POD modes, as shown in figure 5(d). This is expected because POD modes span the optimal linear subspace in the $\mathcal {L}_2$ sense and therefore represent the lower bound for the response projection error. Nevertheless, for this system, the error of the projection onto the EOFs follows this lower bound closely. Importantly, this performance is achieved without any data snapshots and requires only knowledge of the mean flow, as opposed to the long sequence of high-fidelity snapshots needed to obtain converged POD modes. These results are not surprising since resolvent analysis has been very successful to model wall-bounded turbulent flows (McKeon Reference McKeon2017), and as explained in § 2, the EOFs and the combination of all resolvent response modes at all frequencies span the same subspace.

5.2. Sensor placement for reconstruction

We now examine the performance of interpolatory projectors to reconstruct velocity fluctuations in the minimal channel flow from sparse sensor measurements. Again, sensor locations are selected using the pivoted QR approach from Manohar et al. (Reference Manohar, Brunton, Kutz and Brunton2018), and are thus tailored to the specific basis used. The sensor placement distribution reveals the spatial regions that are most informative to reconstruct the dominant coherent structures. Since the state vector of the system contains the three velocity components, sensors are also selecting which component to measure.

Leveraging both of our bases, of POD modes and EOFs, we find that sensors accumulate towards the walls as more of them are added, as shown in figure 6(a). A detailed view of the sensor distribution in the wall-normal direction, presented in figure 6(b), reveals that most of the sensors tailored to the EOFs are positioned closer to the wall and fall within the buffer region. For both bases, the great majority of the sensors measure streamwise velocity, which is expected since that is the most energetic component. Interestingly, there is more variety in the sensors tailored to the POD basis, including more measurements of the spanwise and wall-normal velocities, as shown in figure 6(b).

Figure 6. (a) Spatial locations of sensors designed to reconstruct the velocity fluctuations in minimal channel flow at $Re_{\tau }=185$. Sensors selected via pivoted QR (Manohar et al. Reference Manohar, Brunton, Kutz and Brunton2018) are tailored to the leading $r$ POD modes (top) and the leading $r$ EOFs (bottom). Sensors can measure streamwise (black circles), spanwise (cyan squares) or wall-normal (magenta triangles) velocity components. (b) Mean velocity profile in wall units (bottom) and distribution of sensors in the wall-normal direction for $r=800$ using POD (top) and EOFs (middle). (c) Normalized output projection error as a function of the projector rank $r$ (left) and of the percentage of sensors over the total number of state variables (right). Three interpolatory projectors are considered for each basis, using: $r$ tailored sensors complemented with $2r$ random sensors (dashed), $3r$ random sensors (dash-dotted), and all possible sensors (solid), which is equivalent to an orthogonal projector.

A total of six projectors are built, three for each basis, and the reconstruction $\mathcal {L}_2$ error over the entire sequence of DNS snapshots is evaluated as a function of the projection rank $r$, as shown in figure 6(c). Considering the leading $r$ elements in each basis, we build interpolatory output projectors using $r$ tailored sensors complemented with an additional $2r$ random sensors. Because this is a high-dimensional system with a total of $n\approx 1.6\times 10^5$ states, complementing tailored sensors with random ones is a simple strategy to improve the resulting reconstruction that was suggested in the work of Manohar et al. (Reference Manohar, Brunton, Kutz and Brunton2018). This is implemented by simply taking the leading $2r$ columns of a matrix containing all remaining possible sensors after a random column permutation. These columns are then concatenated to the sampling matrix $\boldsymbol {P}$ used to build the projector. For the sake of comparison, we also build interpolatory projectors using $3r$ random sensors and orthogonal projectors for the subspaces spanned by the POD modes and by the eigenvectors of $\boldsymbol {W_c}$. Orthogonal projectors provide the lower bound in reconstruction error for interpolatory projectors onto the same subspace. Leveraging the EOFs, the reconstruction using tailored sensors explains a large portion of the energy in the flow field, and even with random sensors, $50\,\%$ of the energy is captured from measurements of $1.5\,\%$ of the state variables, as shown in figure 6(c). Furthermore, this is achieved in a data-free manner, requiring only knowledge of the mean flow.