1 Introduction
The laminar flow state is characterized by a lower friction drag than the turbulent one, which implies less energy consumption for many applications, such as transportation means like trains and aircraft. Therefore, control of laminar–turbulent transition is of great interest in many technical areas. The transition scenario depends on a number of parameters, and an overall picture of these different scenarios can be found in Schmid & Henningson (Reference Schmid and Henningson2001). Transition to turbulence in boundary layer flows where free-stream turbulence has an intensity higher than 
${\approx}1\,\%$ occurs rapidly and bypasses the classical scenario triggered by Tollmien–Schlichting (TS) waves, as shown by Arnal & Juillen (Reference Arnal and Juillen1978). When free-stream turbulence is present, a set of low-frequency vortices (Hultgren & Gustavsson Reference Hultgren and Gustavsson1981; Hunt & Durbin Reference Hunt and Durbin1999; Zaki & Saha Reference Zaki and Saha2009; Zhang et al. Reference Zhang, Zaki, Sherwin and Wu2011) enter the boundary layer and causes the appearance of elongated streaky structures of alternating high and low streamwise velocity. The effects of free-stream turbulence on the boundary layer flow were firstly observed in the experimental studies of Klebanoff (Reference Klebanoff1971). The amplitude of such velocity fluctuations grows linearly along the streamwise direction (Andersson, Berggren & Henningson Reference Andersson, Berggren and Henningson1999; Luchini Reference Luchini2000) and is accompanied by growing secondary fluctuations of the streaky structures on the planes perpendicular to the streamwise direction. When the amplitude of such secondary cross-flow fluctuations is sufficiently high turbulent spots appear (Brandt & Henningson Reference Brandt and Henningson2002; Ricco, Luo & Wu Reference Ricco, Luo and Wu2011), which grow and merge further downstream and ultimately lead to a fully turbulent flow. This process was observed both in experiments (Matsubara & Alfredsson Reference Matsubara and Alfredsson2001) and simulations (Brandt, Schlatter & Henningson Reference Brandt, Schlatter and Henningson2004). Thus, the boundary layer can be divided into three zones: (i) an upstream zone where there is high level of receptivity and free-stream turbulence triggers disturbances in the boundary layer, (ii) a middle zone where streaky disturbances grow due to the linear lift-up mechanism and (iii) a downstream zone where the flow nucleates turbulent spots which grow and merge as they propagate downstream until the boundary layer becomes fully turbulent.
The boundary layer flow in the middle zone can often be described with sufficient accuracy by the linearized Navier–Stokes (N–S) equations (Schmid & Henningson Reference Schmid and Henningson2001). The possibility of working with a linear system greatly facilitates the application of flow control techniques. The a priori knowledge of the linear behaviour of the TS waves was exploited in the experiments of Thomas (Reference Thomas1983) and in the simulations of Laurien & Kleiser (Reference Laurien and Kleiser1989) to counteract TS waves and delay transition. Similarly, for bypass transition, Jacobson & Reynolds (Reference Jacobson and Reynolds1998) and Hanson et al. (Reference Hanson, Bade, Belson, Lavoie, Naguib and Rowley2014) exploited the linearity of the dynamical system to show the possibility to damp streaky structures. In those works the a priori knowledge of the system dynamics was used to create ad hoc counter disturbances. Such an ad hoc practice lacks in generality and may require tedious testing, therefore it is appealing to apply optimal control theory to flow control problems. The control theory community has produced many reliable and elegant techniques to tackle linear systems. Among the first successful applications of optimal control theory in fluid mechanics are the works of Joshi, Speyer & Kim (Reference Joshi, Speyer and Kim1997), Bewley & Liu (Reference Bewley and Liu1998), Högberg & Bewley (Reference Högberg and Bewley2000) and Högberg, Bewley & Henningson (Reference Högberg, Bewley and Henningson2003), where optimal control methods are applied to linearized systems and used in fully nonlinear channel flows. More recently, Monokrousos et al. (Reference Monokrousos, Brandt, Schlatter and Henningson2008) showed the successful application of the linear quadratic Gaussian regulator (LQG) for control of streaks triggered by the free-stream turbulence.
In optimal control techniques the final goal is to find the function that takes measurements as input and gives actuation signals as output while minimizing an objective function. Particularly, in classical optimal control methods the optimal solution for linear time-invariant systems is given by solving an algebraic Riccati equation, which consists of a matrix equation whose dimensions are approximately those of the original linear system to control. If the original linear system has large dimensions, as is the case in fluid mechanics, the solution of the algebraic Riccati equation may be extremely computationally demanding. A possible solution is reducing the order of the optimal control problem by keeping only the information useful for the control. This is the idea behind reduced-order models (ROMs). In fact, measurements usually contain only a portion of the total information present in the system and actuators can usually excite only certain structures. In control theory, such limitations posed by sensors and actuators define two properties of the system: its observability and its controllability, respectively. The control problem alone needs only the portion of the system that is observable and controllable. The practice of model reduction in flow control was treated in Bagheri, Brandt & Henningson (Reference Bagheri, Brandt and Henningson2009), Semeraro et al. (Reference Semeraro, Bagheri, Brandt and Henningson2011), Poussot-Vassal & Sipp (Reference Poussot-Vassal and Sipp2015) and Yao & Jaiman (Reference Yao and Jaiman2017). The approach was shown to be successful in the sense that the solution to the control problem was nearly unaffected by the use of a ROM. A classic technique for achieving a ROM is the eigensystem realization algorithm (ERA) (Juang & Pappa Reference Juang and Pappa1985; Semeraro et al. Reference Semeraro, Bagheri, Brandt and Henningson2013a). ERA is based on a set of impulse responses from each input (actuators and disturbances) to each output (measurements). In order to avoid confusion it is noticeable that this description of the ERA is also known as ERA-POD, with POD standing for proper orthogonal decomposition.
Ma, Ahuja & Rowley (Reference Ma, Ahuja and Rowley2011) showed that the ROM achieved by the ERA-POD is equivalent to that achieved by approximate balanced POD truncation. This means that the ROM resulting from the ERA-POD is a projection of the original system onto the set of modes given by the intersection of the set of the most observable flow structures and the set of the most controllable flow structures. Qualitatively, an observable structure is one that generates non-zero outputs whereas a controllable structure is one that can be excited by the inputs. The term ‘most’ is obviously case dependent. For controllable structures it represents the number of flow structures used to recreate with an acceptable small error the flow field obtained by an impulse response. The same reasoning holds for the observable flow structures but with respect to the adjoint system (see Bagheri et al. (Reference Bagheri, Brandt and Henningson2009) for a more detailed discussion of controllability and observability in fluid mechanical systems). Examples of ERA-POD applications in fluid mechanics are found in Semeraro et al. (Reference Semeraro, Bagheri, Brandt and Henningson2013a), for control of a three-dimensional nonlinear TS wave packet, and Sasaki et al. (Reference Sasaki, Morra, Fabbiane, Cavalieri, Hanifi and Henningson2018a), for control of three-dimensional TS waves arising from stochastic disturbances.
A consistent modelling of the inputs implies correctly modelling the space spanned by the disturbances, which may require the use of a basis with as many degrees of freedom as the dimensions of the desired space. Thus, in case the space spanned by a disturbance or an actuator has large dimensions it may become unfeasible to collect all the impulse responses. A similar issue may also happen when the number of outputs is very high. Another possibility to avoid demanding computations for solving the control problem is dropping the use of model-based methods, as discussed by Fabbiane et al. (Reference Fabbiane, Semeraro, Bagheri and Henningson2014), who made use of a learning algorithm that only needs the modelling of the transfer function (TF) from the actuators to the measurements.
The present work addresses the delay of bypass transition in a framework which can be reproduced in wind-tunnel experiments. The disturbance used is stochastic and has high dimensions in order to model the effects of free-stream turbulence in the boundary layer altogether. The number, the location and the type of devices modelled as actuators or sensors are chosen by assuming limitations which occur in experiments. We use a finite number of localized near-wall actuators that resemble ring plasma actuators (Kim & Choi Reference Kim and Choi2016; Kim, Forte & Choi Reference Kim, Forte and Choi2017; Shahriari, Kollert & Hanifi Reference Shahriari, Kollert and Hanifi2018) and localized wall-shear-stress sensors. The techniques to build the ROM and design the controller are chosen by assuming the availability of data which are feasible to collect in experiments. We make use of the ERA-POD, a data-driven algorithm, to generate the ROM for the design of the controller. However, ERA-POD is based on impulse responses from each input to each output, and if the number of inputs is high, as occurs when modelling free-stream turbulence as a disturbance, performing the ERA-POD becomes computationally heavy. Moreover, it is not possible to collect impulse responses from free-stream turbulence in an experimental set-up. So, we present a new method to identify the relevant effects of free-stream turbulence for control purposes and generate a set of impulse responses for the application of the ERA-POD. This new set of impulse responses is smaller than the one needed if the complete disturbance had been used and is based only on measurements. Therefore, the proposed method reduces the computational cost for control design, when compared to the computational cost of the standard application of the ERA-POD in the presence of a large number of inputs, and is feasible in experiments. We design the LQG regulator on the state-space ROM and show that the delay of bypass transition achieved in these realistic conditions is at least as large as the one presented in more idealized studies which do not account for the limitations present in experiments (Monokrousos et al. Reference Monokrousos, Brandt, Schlatter and Henningson2008). Moreover, in order to further discuss the advantages of applying the LQG, which is not feasible in this flow case without the availability of a ROM of reasonable dimensions, we compare the behaviour of the implemented LQG with that of a different optimal control method, the inverse feed-forward control (IFFC), whose action consists of wave cancellation (Sasaki et al. Reference Sasaki, Morra, Fabbiane, Cavalieri, Hanifi and Henningson2018a).
In summary, the contribution of this work consists of: (i) presenting a data-driven method based on localized wall measurements and actuators for the closed-loop control of bypass transition, which exploits the characteristics of the boundary layer flow to model free-stream turbulence in a ROM via ERA-POD and which can be implemented in experiments; (ii) providing evidence via a numerical study that in a realistic framework, considering the limitations which occur in experiments, it is possible to achieve as large delay of bypass transition as that obtained in more idealized cases found in the literature.
The paper is structured as follows. In § 2 the equations used to describe the full or reduced system dynamics are introduced; in § 3 the control techniques of interest are briefly described; in § 4 the details of the framework for the nonlinear simulations are outlined; in § 5 the identification techniques and the used identified models are presented; in § 6 the behaviour of the designed controller in the nonlinear N–S equations is assessed. A summary of the main conclusions is given in § 7.
2 Governing equations
2.1 Dynamical system
The N–S equations can be expressed in terms of the perturbation quantities as
$$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{q}^{\prime }}{\unicode[STIX]{x2202}t}}=-(\boldsymbol{q}^{\prime }\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{q}_{B}-(\boldsymbol{q}_{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{q}^{\prime }-(\boldsymbol{q}^{\prime }\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{q}^{\prime }-\unicode[STIX]{x1D735}p^{\prime }+Re^{-1}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{q}^{\prime }+\boldsymbol{f}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{q}^{\prime }=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{q}^{\prime }=\boldsymbol{q}_{0}^{\prime }\quad \text{at }t=0, & \displaystyle\end{eqnarray}$$
$\boldsymbol{q}^{\prime }=\boldsymbol{q}^{\prime }(\boldsymbol{x},t)$ is the perturbation velocity vector, 
$\boldsymbol{q}_{B}=\boldsymbol{q}_{B}(\boldsymbol{x})$ the unperturbed velocity vector, 
$p^{\prime }=p^{\prime }(\boldsymbol{x},t)$ the perturbation pressure, 
$\boldsymbol{f}=\boldsymbol{f}(\boldsymbol{x},t)$ a body force vector, 
$Re$ the Reynolds number, 
$t$ the time, 
$\boldsymbol{x}=(x_{1},x_{2},x_{3})^{T}$ the space variable and 
$\unicode[STIX]{x1D735}=(\unicode[STIX]{x2202}_{x_{1}},\unicode[STIX]{x2202}_{x_{2}},\unicode[STIX]{x2202}_{x_{3}})$ the gradient operator. Here, the unperturbed velocity vector 
$\boldsymbol{q}_{B}$ is the solution of an evolving zero-pressure-gradient Blasius boundary layer. The velocity perturbation 
$\boldsymbol{q}^{\prime }$ satisfies no-slip conditions at the wall 
$x_{2}=0$ and Neumann conditions at free stream 
$x_{2}=L_{x_{2}}$. Periodicity is assumed along the spanwise direction 
$x_{3}$ and enforced along the streamwise direction 
$x_{1}$ by means of an artificial forcing 
$\boldsymbol{f}_{BC}=\unicode[STIX]{x1D706}(x_{1})(\boldsymbol{q}_{B}-\boldsymbol{q}^{\prime })$, which is placed in a fringe region at the outlet; 
$\unicode[STIX]{x1D706}(\boldsymbol{x})$ is a non-negative function which is non-zero only within the fringe region; 
$\boldsymbol{f}_{BC}$ forces all perturbations to zero and modifies 
$\boldsymbol{q}_{B}$ to be periodic (see Nordström, Nordin & Henningson (Reference Nordström, Nordin and Henningson1999) for details about 
$\boldsymbol{f}_{BC}$).
The flow control problem consists of finding the correct external action that modifies the fluid dynamics to achieve a specific goal. In our case, such external action can take the form of a boundary condition or a body force and can be expressed as a function of time and space. It follows that the problem can be split into finding the correct spatial distribution of such an action and its time modulation. In the present work it is assumed that the spatial distribution and the time modulation of the external action are decoupled. The spatial distribution is prescribed, so the flow control problem reduces to the computation of its time-varying amplitude. From now on this time-varying scalar is referred to as the input. Using a finite number of actuators 
$N_{u}$, the external action used for control reads
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{f}_{u}=\mathop{\sum }_{k=1}^{N_{u}}u(t)_{k}\boldsymbol{b}(\boldsymbol{x})_{k}, & \displaystyle\end{eqnarray}$$ where 
$\boldsymbol{b}(\boldsymbol{x})_{k}$ is the spatial shape of the 
$k$th body force, and 
$u(t)_{k}$ the corresponding time variation. The latter represents the control input.
Free-stream turbulence is modelled as a forcing in the fringe region; 
$\boldsymbol{f}_{BC}$ is modified to force 
$\boldsymbol{q}^{\prime }$ to be equal to a prescribed perturbation that mimics the presence of free-stream turbulence. The prescribed perturbation is of the form
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{q}_{FST}^{\prime }=\mathop{\sum }_{\unicode[STIX]{x1D6FC}}\mathop{\sum }_{\unicode[STIX]{x1D6FD}}\mathop{\sum }_{\unicode[STIX]{x1D714}}\unicode[STIX]{x1D6F7}(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714})\hat{\boldsymbol{q}}^{\prime }(x_{2},\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714})\text{e}^{\text{i}(\unicode[STIX]{x1D6FC}x_{1}+\unicode[STIX]{x1D6FD}x_{3}-\unicode[STIX]{x1D714}t)}, & \displaystyle\end{eqnarray}$$ with 
$\hat{\boldsymbol{q}}^{\prime }$ an eigensolution to the Orr–Sommerfeld–Squire eigenvalue problem for a parallel flow in a semi-bounded domain, 
$\unicode[STIX]{x1D6FC}$ the streamwise wavenumber, 
$\unicode[STIX]{x1D6FD}$ the spanwise wavenumber and 
$\unicode[STIX]{x1D714}$ the angular frequency. Free-stream disturbances are thus expanded as a sum of eigenfunctions of the linearized parallel-flow problem (see Brandt et al. (Reference Brandt, Schlatter and Henningson2004) for more details).
A linearized version of the N–S equations about 
$\boldsymbol{q}_{B}$ can be obtained by dropping the nonlinear term 
$(\boldsymbol{q}^{\prime }\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{q}^{\prime }$ from (2.1a).
Plant. Computational box (cut along 
$x_{1}$): frame of reference 
$x_{1}x_{2}x_{3}$, sensors 
$y$ and 
$z$ (black circles), actuators 
$u$ (white circles) and controller 
$G^{yu}$. Contour plots: perturbation part of the streamwise velocity, 
$q_{1}^{\prime }$, snapshot of an uncontrolled case at time 
$t=t^{\ast }$; the isolines of the contour plots are all for the same set of values. Boundary layer, 
$\unicode[STIX]{x1D6FF}_{99}$, shown on the left wall of the box.
2.2 Reduced-order dynamical system
The linear dynamical system used for the application of control theory techniques is a ROM and reads
$$\begin{eqnarray}\displaystyle & \displaystyle \dot{\boldsymbol{q}}(t)=\unicode[STIX]{x1D63C}\boldsymbol{q}(t)+\unicode[STIX]{x1D63D}\boldsymbol{u}(t)+\unicode[STIX]{x1D648}_{d}\boldsymbol{d}(t), & \displaystyle\end{eqnarray}$$ where 
$\boldsymbol{q}=\boldsymbol{q}(t)$ is the 
$N\times 1$ state vector (which generally is not exactly the same quantity represented by 
$\boldsymbol{q}^{\prime }$), 
$\dot{\boldsymbol{q}}$ is its time derivative, 
$\unicode[STIX]{x1D63C}$ the 
$N\times N$ matrix that defines the system dynamics, 
$\unicode[STIX]{x1D63D}$ the 
$N\times N_{u}$ matrix that characterizes the control inputs, 
$\boldsymbol{u}=\boldsymbol{u}(t)$ a 
$N_{u}\times 1$ column vector containing all the input amplitudes 
$u(t)_{k}$, 
$\unicode[STIX]{x1D648}_{d}$ the 
$N\times N_{d}$ matrix that characterizes the disturbance inputs and 
$\boldsymbol{d}=\boldsymbol{d}(t)$ a 
$N_{d}\times 1$ column vector containing all the input amplitudes 
$d(t)_{k}$. Here, 
$N$ is the degree of freedom of the ROM, 
$N_{u}$ the number of control inputs and 
$N_{d}$ the number of disturbance inputs.
We also assume to have access to two finite sets of measurements: 
$\boldsymbol{y}(t)$, 
$N_{y}\times 1$ and 
$\boldsymbol{z}(t)$, 
$N_{z}\times 1$, where 
$N_{y}$ and 
$N_{z}$ represent the respective number of measurements available. It holds that
$$\begin{eqnarray}\displaystyle \boldsymbol{y}(t)=\unicode[STIX]{x1D63E}_{y}\boldsymbol{q}(t),\quad \boldsymbol{z}(t)=\unicode[STIX]{x1D63E}_{z}\boldsymbol{q}(t), & & \displaystyle\end{eqnarray}$$ where the 
$N_{y}\times N$ matrix 
$\unicode[STIX]{x1D63E}_{y}$ and the 
$N_{z}\times N$ matrix 
$\unicode[STIX]{x1D63E}_{z}$ characterize the measurements in the ROM. From now on 
$\boldsymbol{y}(t)$ and 
$\boldsymbol{z}(t)$ are referred to as outputs. Equations (2.4) and (2.5) form a ROM state-space representation of the system.
A different description of the system can be given by means of transfer functions. TFs are built by performing the Laplace transform on the state-space representation and in general describe the system as a function of the angular frequency 
$\unicode[STIX]{x1D714}$ only. Here, sensors and actuators are placed on straight lines along the spanwise direction (figure 1), with 
$N_{u}=N_{y}=N_{z}$. Since the flow is periodic in the spanwise direction, TFs, inputs and outputs can be expressed as functions of the spanwise wavenumber 
$\unicode[STIX]{x1D6FD}$ as well. The description by means of TFs reads
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}{\hat{y}}(\unicode[STIX]{x1D714},\unicode[STIX]{x1D6FD}_{k})={\hat{G}}^{uy}(\unicode[STIX]{x1D714},\unicode[STIX]{x1D6FD}_{k})\hat{u} (\unicode[STIX]{x1D714},\unicode[STIX]{x1D6FD}_{k})+{\hat{G}}^{dy}(\unicode[STIX]{x1D714},\unicode[STIX]{x1D6FD}_{k})\hat{d}(\unicode[STIX]{x1D714},\unicode[STIX]{x1D6FD}_{k}),\\ \hat{z}(\unicode[STIX]{x1D714},\unicode[STIX]{x1D6FD}_{k})={\hat{G}}^{uz}(\unicode[STIX]{x1D714},\unicode[STIX]{x1D6FD}_{k})\hat{u} (\unicode[STIX]{x1D714},\unicode[STIX]{x1D6FD}_{k})+{\hat{G}}^{dz}(\unicode[STIX]{x1D714},\unicode[STIX]{x1D6FD}_{k})\hat{d}(\unicode[STIX]{x1D714},\unicode[STIX]{x1D6FD}_{k}),\end{array}\right\} & & \displaystyle\end{eqnarray}$$ where 
${\hat{G}}={\hat{G}}(\unicode[STIX]{x1D714},\unicode[STIX]{x1D6FD}_{k})$ are the TFs, 
${\hat{y}}={\hat{y}}(\unicode[STIX]{x1D714},\unicode[STIX]{x1D6FD}_{k})$ and 
$\hat{z}=\hat{z}(\unicode[STIX]{x1D714},\unicode[STIX]{x1D6FD}_{k})$ the outputs in the frequency domain, 
$\hat{u} =\hat{u} (\unicode[STIX]{x1D714},\unicode[STIX]{x1D6FD}_{k})$ and 
$\hat{d}=\hat{d}(\unicode[STIX]{x1D714},\unicode[STIX]{x1D6FD}_{k})$ the inputs in the frequency domain and 
$k$ is used to stress the fact that the number of outputs is finite, so there is a finite number of available wavenumbers. From now on all the variables denoted by a hat symbol are function of 
$(\unicode[STIX]{x1D714},\unicode[STIX]{x1D6FD}_{k})$, and the explicit writing 
$(\unicode[STIX]{x1D714},\unicode[STIX]{x1D6FD}_{k})$ is dropped.
The description of the system by means of TFs can be translated in the physical domain by performing the inverse Fourier transform to (2.6), leading to
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle y(t)_{k}=\int _{0}^{t}\mathop{\sum }_{m=1}^{N_{u}}G_{km}^{uy}(t-\unicode[STIX]{x1D70F})u(\unicode[STIX]{x1D70F})_{m}\,\text{d}\unicode[STIX]{x1D70F}+\int _{0}^{t}\mathop{\sum }_{m=1}^{N_{d}}G_{km}^{dy}(t-\unicode[STIX]{x1D70F})d(\unicode[STIX]{x1D70F})_{m}\,\text{d}\unicode[STIX]{x1D70F},\\[18.0pt] \displaystyle z(t)_{k}=\int _{0}^{t}\mathop{\sum }_{m=1}^{N_{u}}G_{km}^{uz}(t-\unicode[STIX]{x1D70F})u(\unicode[STIX]{x1D70F})_{m}\,\text{d}\unicode[STIX]{x1D70F}+\int _{0}^{t}\mathop{\sum }_{m=1}^{N_{d}}G_{km}^{dz}(t-\unicode[STIX]{x1D70F})d(\unicode[STIX]{x1D70F})_{m}\,\text{d}\unicode[STIX]{x1D70F},\end{array}\right\} & & \displaystyle\end{eqnarray}$$ with 
$k$ the output index, and 
$m$ the input index. Equation (2.7) can also be obtained by substituting the solution to (2.4), with 
$\boldsymbol{q}(0)=0$, into (2.5). This gives the identities
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}G_{km}^{uy}(t-\unicode[STIX]{x1D70F})=\unicode[STIX]{x1D63E}_{y,k}\text{e}^{\unicode[STIX]{x1D63C}(t-\unicode[STIX]{x1D70F})}\unicode[STIX]{x1D63D}_{m},\quad G_{km}^{dy}(t-\unicode[STIX]{x1D70F})=\unicode[STIX]{x1D63E}_{y,k}\text{e}^{\unicode[STIX]{x1D63C}(t-\unicode[STIX]{x1D70F})}\unicode[STIX]{x1D648}_{d,m},\\[6.0pt] G_{km}^{uz}(t-\unicode[STIX]{x1D70F})=\unicode[STIX]{x1D63E}_{z,k}\text{e}^{\unicode[STIX]{x1D63C}(t-\unicode[STIX]{x1D70F})}\unicode[STIX]{x1D63D}_{m},\quad G_{km}^{dz}(t-\unicode[STIX]{x1D70F})=\unicode[STIX]{x1D63E}_{z,k}\text{e}^{\unicode[STIX]{x1D63C}(t-\unicode[STIX]{x1D70F})}\unicode[STIX]{x1D648}_{d,m},\end{array}\right\} & & \displaystyle\end{eqnarray}$$ where 
$\unicode[STIX]{x1D63E}_{y,k}$ and 
$\unicode[STIX]{x1D63E}_{z,k}$ the 
$k$th row of 
$\unicode[STIX]{x1D63E}_{y}$ and 
$\unicode[STIX]{x1D63E}_{z}$, and 
$\unicode[STIX]{x1D63D}_{m}$ and 
$\unicode[STIX]{x1D648}_{d,m}$ the 
$m$th column of 
$\unicode[STIX]{x1D63D}$ and 
$\unicode[STIX]{x1D648}_{d}$.
3 Control techniques
The present configuration of outputs and inputs together with the convective nature of the flow make all the control techniques described in this section feed-forward configurations (Belson et al. Reference Belson, Semeraro, Rowley and Henningson2013).
The control techniques used in the present work are all based on the assumption that the input 
$u(t)_{k}$ is a function of the upstream outputs 
$y(t)_{m}$,
$$\begin{eqnarray}\displaystyle & \displaystyle u(t)_{k}=\int _{0}^{t}\mathop{\sum }_{m=1}^{N_{y}}G_{km}^{yu}(t-\unicode[STIX]{x1D70F})y(\unicode[STIX]{x1D70F})_{m}\,\text{d}\unicode[STIX]{x1D70F}, & \displaystyle\end{eqnarray}$$ with 
$k=1,2,\ldots ,N_{u}$; 
$G_{km}^{yu}$ in (3.1) differs from the quantities in (2.7) because it does not describe the open-loop input–output dynamics. It is designed via control theory for the prescribed closed-loop dynamics. Since 
$\boldsymbol{q}_{B}$ is independent of the spanwise direction 
$x_{3}$, the instantaneous linearized system dynamics is homogeneous along 
$x_{3}$. The latter and the fact that the outputs are all given by the same type of sensor allows us to drop the usage of the index 
$k$ in (3.1) to have 
$G_{m}^{yu}$. Then, (3.1) can be rewritten as
$$\begin{eqnarray}\displaystyle & \displaystyle u(t)_{k}=\int _{0}^{t}\mathop{\sum }_{m=1}^{N_{y}}G_{m}^{yu}(t-\unicode[STIX]{x1D70F})y(\unicode[STIX]{x1D70F})_{m+k-1}\,\text{d}\unicode[STIX]{x1D70F}, & \displaystyle\end{eqnarray}$$ where for 
$m+k-1>N_{y}$ spanwise periodicity implies the use of 
$m+k-1-N_{y}$.
3.1 Linear quadratic Gaussian regulator
The technique is based on a linear model, aims at minimizing a quadratic cost function and assumes the presence of Gaussian white noise disturbances.
Gaussian white noise is added on the output 
$\boldsymbol{y}(t)$ in the ROM, which reads 
$\boldsymbol{y}(t)=\unicode[STIX]{x1D63E}_{y}\boldsymbol{q}(t)+\boldsymbol{n}(t)$, with 
$\boldsymbol{n}(t)$, 
$N_{y}\times 1$ the time modulation of the noise. Here, 
$\boldsymbol{d}(t)$ in (2.4) is also treated as white noise. There is no addition of noise to the output 
$\boldsymbol{z}(t)$ because it represents a reference output to minimize, whose measurement is not available in reality. The noise on the output 
$\boldsymbol{y}(t)$ corresponds to noise in a real available measurement.
The covariance matrices associated with 
$\boldsymbol{d}(t)$ and 
$\boldsymbol{n}(t)$ are 
$\unicode[STIX]{x1D651}_{d}$, 
$N_{d}\times N_{d}$, and 
$\unicode[STIX]{x1D651}_{n}$, 
$N_{y}\times N_{y}$, respectively, and are both diagonal and constant because of the assumption that 
$\boldsymbol{d}(t)$ and 
$\boldsymbol{n}(t)$ are white noise disturbances; in particular, the following can be written
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D651}_{d}=v_{d}\unicode[STIX]{x1D644},\quad \unicode[STIX]{x1D651}_{n}=v_{n}\unicode[STIX]{x1D644}, & & \displaystyle\end{eqnarray}$$ with 
$v_{d}>0$ and 
$v_{n}>0$ real scalars and 
$\unicode[STIX]{x1D644}$ the identity matrix.
The technique consists of finding 
$G_{m}^{yu}$ by minimizing a prescribed 
${\mathcal{H}}_{2}$-norm of interest. The disturbances 
$\boldsymbol{d}(t)$ and 
$\boldsymbol{n}(t)$ are treated as random variables, so the objective function of interest is defined as the expected value of an 
${\mathcal{H}}_{2}$-norm. Here, the objective function contains both the reference output 
$\boldsymbol{z}(t)$ and the input for the control 
$\boldsymbol{u}(t)$, which is added to avoid an infinite amplitude of the input signal, penalizing excessive control action. The objective function reads
$$\begin{eqnarray}\displaystyle & \displaystyle J=\mathbb{E}\left[\lim _{T\rightarrow \infty }\frac{1}{T}\int _{0}^{T}\boldsymbol{z}(t)^{T}\unicode[STIX]{x1D64C}\boldsymbol{z}(t)+\boldsymbol{u}(t)^{T}\unicode[STIX]{x1D64D}\boldsymbol{u}(t)\,\text{d}t\right], & \displaystyle\end{eqnarray}$$ where the 
$N\times N$ matrices 
$\unicode[STIX]{x1D64C}$ and 
$\unicode[STIX]{x1D64D}$ are design weights. The operator 
$\mathbb{E}[\bullet ]$ represents the expected value. From now on the matrices 
$\unicode[STIX]{x1D651}_{d}$, 
$\unicode[STIX]{x1D651}_{n}$, 
$\unicode[STIX]{x1D64C}$ and 
$\unicode[STIX]{x1D64D}$ are referred to as weight matrices or design weights.
In the LQG it is assumed that 
$\boldsymbol{u}(t)$ is a linear function of the states, but it is also assumed that not all the states are known at each time instant, so a second system for state estimation is introduced. The estimation system makes use of the known outputs to reconstruct the states at each instant of time, and is designed to minimize the estimation error. Thus, in addition to the minimization of the objective function (3.4) to compute the input that controls the system, the estimation introduces a second minimization problem. Generally these two minimization problems are coupled, but in the LQG they are independent and solved separately. They consist of the linear quadratic regulator, which solves the control problem by assuming full-state information, and the Kalman filter, which solves the estimation problem by assuming stochastic disturbances on the outputs. The solution of the LQG is the combination of the two independent solutions. More details about the LQG are given in appendix A, whereas a thorough description can be found in Lewis & Syrmos (Reference Lewis and Syrmos1995).
3.2 Inversion feed-forward control
Inversion feed-forward control is a technique developed in the frequency domain, and is based on a system described by TFs (2.6). The technique is exactly the same one used in Sasaki et al. (Reference Sasaki, Morra, Fabbiane, Cavalieri, Hanifi and Henningson2018a), but in that work it is referred to as wave cancellation. The authors decided to change the nomenclature to adopt the name used in the control community.
The contribution of the disturbance 
$\hat{d}$ in the second equation of (2.6) may also be expressed as
$$\begin{eqnarray}\displaystyle & \displaystyle {\hat{G}}^{dz}\hat{d}={\hat{G}}^{yz}{\hat{y}}+\hat{p}, & \displaystyle\end{eqnarray}$$ where 
${\hat{G}}^{yz}$ is a TF to design in order to maximize the extraction of information from 
${\hat{y}}$, while 
$\hat{p}$ is the residual part of the information in 
$\hat{z}$ which is not retrieved by 
${\hat{G}}^{yz}{\hat{y}}$. The loss of information, i.e. 
$\hat{p}\neq 0$, may be unavoidable and can be seen by resorting to the state-space representation. The matrix 
$\unicode[STIX]{x1D63E}_{y}$ that characterizes the output 
$\boldsymbol{y}(t)$ does not necessarily span the same space spanned by the matrix that describes the system dynamics 
$\unicode[STIX]{x1D63C}$, so the outputs 
$\boldsymbol{y}(t)$, being in a subspace, cannot reconstruct the whole state space. It follows that the only portion of the signal 
$\hat{z}(\unicode[STIX]{x1D714},\unicode[STIX]{x1D6FD}_{k})$ which can be obtained from the outputs 
$\boldsymbol{y}(t)$ is
$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\hat{z}}={\hat{G}}^{uz}\hat{u} +{\hat{G}}^{yz}{\hat{y}}, & \displaystyle\end{eqnarray}$$ where 
$\tilde{\hat{z}}$ is an estimate of 
$\hat{z}$.
The objective of the control problem is the annihilation of the output 
$\tilde{\hat{z}}$. Then, a straightforward strategy to solve the problem is to impose 
$\tilde{\hat{z}}=0$, which is the basic idea behind IFFC. Assuming 
$\hat{u} =\hat{K}{\hat{y}}$ in (3.6) gives
$$\begin{eqnarray}\displaystyle \hat{K}=({\hat{G}}^{uz})^{-1}{\hat{G}}^{yz}, & & \displaystyle\end{eqnarray}$$ where 
$\hat{K}$ solves the control problem in the frequency–wavenumber domain. The result in (3.7) is ill conditioned in the zeros of 
${\hat{G}}^{uz}$, which may lead to spurious high amplitudes of the input. Moreover, model uncertainties are not considered, and unstable zeros in 
${\hat{G}}^{uz}$ would lead to an unstable controller, which is rarely appreciated in practice. Such limitations are addressed in Devasia (Reference Devasia2002), where the TFs, the inputs and the outputs are functions of the angular frequency 
$\unicode[STIX]{x1D714}$ only. Here, the same approach is used with some modification to account for a system description as a function of 
$\unicode[STIX]{x1D714}$ and 
$\unicode[STIX]{x1D6FD}_{k}$, following Sasaki et al. (Reference Sasaki, Morra, Fabbiane, Cavalieri, Hanifi and Henningson2018a). The technique makes use of two weights, 
$\hat{R}$ and 
$\hat{Q}$, which here are taken as constant, and solves the control problem by minimizing the following prescribed objective function
$$\begin{eqnarray}\displaystyle & \displaystyle J=\int _{-\infty }^{\infty }\mathop{\sum }_{k}\left(\hat{u} ^{H}\hat{R}\hat{u} +\hat{z}^{H}\hat{Q}\hat{z}\right)\unicode[STIX]{x0394}\unicode[STIX]{x1D6FD}_{k}\,\text{d}\unicode[STIX]{x1D714}, & \displaystyle\end{eqnarray}$$ where the superscript 
$H$ indicates the complex conjugate transpose. The presence of the objective function turns the nature of the problem into an 
${\mathcal{H}}_{2}$ optimal control problem, whose solution is given by
$$\begin{eqnarray}\displaystyle & \displaystyle \hat{K}=\frac{({\hat{G}}^{uz})^{H}\hat{Q}{\hat{G}}^{yz}}{\hat{R}+({\hat{G}}^{uz})^{H}\hat{Q}{\hat{G}}^{uz}}. & \displaystyle\end{eqnarray}$$ The inverse Fourier transform of 
$\hat{K}$ gives 
$G_{m}^{yu}$ as in (3.2); only the causal part of 
$G_{m}^{yu}$ is used for control, since actuation must be decided based solely on present or past information from the sensors.
This technique was already applied by Sasaki et al. (Reference Sasaki, Morra, Fabbiane, Cavalieri, Hanifi and Henningson2018a,Reference Sasaki, Tissot, Cavalieri, Silvestre, Jordan and Biaub) and shown to be successful in the control of Kelvin–Helmholtz and Tollmien–Schlichting waves, where the equivalence between LQG and IFFC for the damping of TS waves was also shown.
4 Plant
The domain of interest is a box as shown in figure 1, where the white symbols represent the outputs and the black symbols the inputs. For flow simulations the pseudo-spectral code SIMSON (Chevalier et al. Reference Chevalier, Schlatter, Lundbladh and Henningson2007) is used. Here, the reference length is taken to be the displacement thickness of the boundary layer at the inlet 
$\unicode[STIX]{x1D6FF}_{0}^{\ast }$ and the reference velocity is the free-stream velocity 
$U_{\infty }$. In all of the present simulations the Reynolds number is 
$Re=U_{\infty }\unicode[STIX]{x1D6FF}_{0}^{\ast }/\unicode[STIX]{x1D708}=300$. All the results that involve transition to turbulence are performed by means of large-eddy simulation (LES) on a box of dimensions 
$(L_{x_{1}},L_{x_{2}},L_{x_{3}})=(4000,60,50)$ with 
$(N_{x_{1}},N_{x_{2}},N_{x_{3}})=(1024,121,108)$ points for the discretization. The effect of the LES filter (see Schlatter, Stolz & Kleiser (Reference Schlatter, Stolz and Kleiser2004), Schlatter, Stolz & Kleiser (Reference Schlatter, Stolz, Kleiser, Lamballais, Friedrich, Geurts and Métais2006a) and Schlatter, Stolz & Kleiser (Reference Schlatter, Stolz and Kleiser2006b) for details) in the area where the flow dynamics is linear is negligible (Monokrousos et al. Reference Monokrousos, Brandt, Schlatter and Henningson2008). All the results that do not need to include the fully turbulent regime are performed by direct-numerical simulations (DNS) on a box of dimensions 
$(L_{x_{1}},L_{x_{2}},L_{x_{3}})=(1000,60,50)$ with 
$(N_{x_{1}},N_{x_{2}},N_{x_{3}})=(1152,121,108)$ points for the discretization. The points along the wall-parallel directions are equi-spaced, whereas along the wall normal there are Gauss–Lobatto points.
The free-stream turbulence is modelled by superposition of 200 random modes from the continuous Orr–Sommerfeld–Squire spectrum. The integral length scale and the turbulent intensity used for all presented results are respectively 
$L=7.5\unicode[STIX]{x1D6FF}_{0}^{\ast }$ and 
$Tu=3.0\,\%$, considering the free-stream turbulence spectrum in Brandt et al. (Reference Brandt, Schlatter and Henningson2004).
As shown in figure 1, the input and output devices are placed along straight lines. The first set of outputs is placed at 
$x_{1,y}=250$, which is downstream of the zone with high receptivity. The second set of outputs is placed at 
$x_{1,z}=400$, since after that position the nonlinearities start to be non-negligible. Input signals, corresponding to actuators, are generated at 
$x_{1,u}=325$ to have the same 
$\unicode[STIX]{x0394}x_{1}$ between input and outputs, such that the travelling time of a disturbance from the first set of outputs to the inputs is roughly the same as the travelling time from the inputs to the second set of outputs. The chosen location for the devices is also optimal in terms of identification accuracy for control design, as shown in appendix B. The number of devices along the spanwise direction is the same for each set and it is equal to 
$N_{u}=N_{y}=N_{z}=36$. Such a choice is motivated by analysing the wavenumber spectrum of the average disturbance energy. According to Shannon information theorem the sampling wavenumber needs to be at least twice the wavenumber of interest. In this case measuring the highest non-negligible spanwise fluctuation would require at least 18 devices. In order to have a better measurement of the spanwise fluctuations 36 devices are used. The devices are equi-spaced along the spanwise direction. The shape of the input actuator is given by
$$\begin{eqnarray}\displaystyle \boldsymbol{b}(\boldsymbol{x})=\{0,b_{2}(\boldsymbol{x}),0\}^{\text{T}} & & \displaystyle\end{eqnarray}$$with
$$\begin{eqnarray}\displaystyle & \displaystyle b_{2}(\boldsymbol{x})=\exp \left[-{\displaystyle \frac{(x_{1}-x_{1,0})^{2}}{\unicode[STIX]{x1D70E}_{x_{1}}^{2}}}-{\displaystyle \frac{(x_{2})^{2}}{\unicode[STIX]{x1D70E}_{x_{2}}^{2}}}-{\displaystyle \frac{(x_{3}-x_{3,0})^{2}}{\unicode[STIX]{x1D70E}_{x_{3}}^{2}}}\right], & \displaystyle\end{eqnarray}$$ where 
$\unicode[STIX]{x1D70E}_{x_{1}}=3$, 
$\unicode[STIX]{x1D70E}_{x_{2}}=5$ and 
$\unicode[STIX]{x1D70E}_{x_{3}}=1.5$. The actuator shape resembles that of ring plasma actuators (see Kim & Choi Reference Kim and Choi2016; Kim et al. Reference Kim, Forte and Choi2017; Shahriari et al. Reference Shahriari, Kollert and Hanifi2018), generating a body force in the wall-normal direction. This is efficient to excite or cancel streaks due to the lift-up effect. A detailed analysis on the effect of the actuator shape is presented in Sasaki et al. (Reference Sasaki, Morra, Cavalieri, Hanifi and Henningson2019), which is the parallel work to the present one.
The outputs are computed as
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{1}{S}\int _{S}\frac{\unicode[STIX]{x2202}q_{1}^{\prime }}{\unicode[STIX]{x2202}x_{2}}\bigg|_{x_{2}=0}\,\text{d}S, & \displaystyle\end{eqnarray}$$ with 
$S$ the area on the wall where the measure is taken. This is an averaged measure of the shear stress associated with the perturbation part of the streamwise velocity component on the wall.
5 Reduced-order modelling and control design
Both control methods introduced in § 3 are model based. The IFFC technique requires knowledge of two TFs, 
${\hat{G}}^{zu}$ and 
${\hat{G}}^{yz}$; 
${\hat{G}}^{zu}$ is by definition the Fourier transform of the output signal resulting from an impulse-response simulation of the linearized N–S equations, whereas 
${\hat{G}}^{yz}$ needs to be modelled. The LQG technique, instead, requires the knowledge of the matrices 
$\unicode[STIX]{x1D63C},~\unicode[STIX]{x1D63D},~\unicode[STIX]{x1D63E}_{y},~\unicode[STIX]{x1D63E}_{z}$ and 
$\unicode[STIX]{x1D648}_{d}$, which characterize the ROM and need to be modelled.
The techniques used for this modelling are introduced in the remainder of this section, and are all based on input–output data, which are usually available in experiments. In input–output data part of the information about the system dynamics is lost. However, its usage is a reasonable design choice, since the control techniques work only with the observable and controllable structures, whose time evolution is described by input–output signals. In fact, by definition, the information lost in input–output data is that associated with the unobservable and uncontrollable structures.
The present configuration of outputs and inputs together with the convective nature of the flow allow us to estimate the downstream outputs 
$\boldsymbol{z}(t)$ from the upstream outputs 
$\boldsymbol{y}(t)$. This fact is exploited in the following part of this section.
5.1 Empirical TFs
The estimation of downstream outputs 
$\hat{z}$ by means of upstream outputs 
${\hat{y}}$ can be performed by designing a TF 
${\hat{G}}^{yz}$. Here, 
${\hat{G}}^{yz}$ is computed by means of an identification technique using the information extracted from the output data. The TF obtained in this way is referred to as empirical TF. This method was introduced in Sasaki et al. (Reference Sasaki, Piantanida, Cavalieri and Jordan2017) for the estimation of a turbulent jet and applied in Sasaki et al. (Reference Sasaki, Tissot, Cavalieri, Silvestre, Jordan and Biau2018b) for the closed-loop control of a two-dimensional shear layer. Here, the approach is extended to a flow with spanwise periodicity, i.e. outputs are functions of 
$\unicode[STIX]{x1D6FD}_{k}$ as well. It was shown in Bendat & Piersol (Reference Bendat and Piersol2011) that the optimal frequency response, in the least square sense, is defined from the auto- and cross-spectra of the input and output signals
$$\begin{eqnarray}\displaystyle & \displaystyle {\hat{G}}_{yz}=\frac{{\hat{S}}_{yz}}{{\hat{S}}_{yy}}, & \displaystyle\end{eqnarray}$$ where 
${\hat{S}}_{yy}$ and 
${\hat{S}}_{yz}$ are respectively the auto- and cross-spectra of the input and output signals. Both 
${\hat{S}}_{yy}$ and 
${\hat{S}}_{yz}$ are computed as the expected values of 
${\hat{y}}^{H}{\hat{y}}$ and 
${\hat{y}}^{H}\hat{z}$, which are obtained via the process of ensemble averaging (Bendat & Piersol Reference Bendat and Piersol2011). Equation (5.1), sometimes referred to as an 
${\hat{H}}_{1}$ estimator (Rocklin, Crowley & Vold Reference Rocklin, Crowley and Vold1985), minimizes the error due to noise in the output.
One desirable property of an 
$H_{1}$ estimator is that the prediction error is linearly uncorrelated to the available output signal (Rocklin et al. Reference Rocklin, Crowley and Vold1985; Bendat & Piersol Reference Bendat and Piersol2011). Any remaining errors correlated to the available signal are either due to the presence of noise in the measurements or to spectral leakage, which is unavoidable because the signal is not exactly periodic in time. Spectral leakage can be minimized by using long time series, by windowing the signal for the ensemble averaging or via the calculation of an improved frequency response, as outlined in the following section.
5.2 Improved frequency response
The method considered here is referred to as improved frequency response and consists of improving the accuracy of a TF by means of an iterative algorithm. This allows us to obtain a more accurate linear approximation of a system and is particularly interesting when the impulse responses of the disturbances are not available or unfeasible to collect, as is the case for the free-stream turbulence or experimental implementations. The method is designed to minimize noise, spectral leakage and capture some nonlinearity (Schoukens, Rolain & Pintelon Reference Schoukens, Rolain and Pintelon1998).
The algorithm is initialized with a first-guess TF 
${\hat{G}}_{0}^{yz}$, which may be, for instance, the result obtained from (5.1). Then, the estimation error, which is the difference between the signal obtained by using 
${\hat{G}}_{0}^{yz}$ and the available output 
${\hat{y}}$, is computed. The error reads
$$\begin{eqnarray}\displaystyle & \displaystyle e(t)_{k}=z(t)_{k}-\int _{0}^{t}\mathop{\sum }_{m}G_{0,m}^{yz}(t-\unicode[STIX]{x1D70F})y(t)_{m}\,\text{d}\unicode[STIX]{x1D70F} & \displaystyle\end{eqnarray}$$ with 
$G_{0,m}^{yz}(t)$ the inverse Fourier transform of 
${\hat{G}}_{0}^{yz}$. Then, the TF between the error and the available output 
${\hat{y}}$ is computed as 
${\hat{G}}_{e}^{yz}={\hat{S}}_{ye}/{\hat{S}}_{yy}$, which is used to update the initial TF as 
${\hat{G}}_{1}^{yz}={\hat{G}}_{0}^{yz}+{\hat{G}}_{e}^{yz}$. Iterations are performed until the error TF is minimized.
5.3 Eigensystem realization algorithm using TFs
In § 3 it was shown that in order to design the LQG regulator it is necessary to solve two algebraic Riccati equations, and the computational power required for their solution grows quickly with the dimensions of the matrix 
$\unicode[STIX]{x1D63C}$, which is the linear time-invariant operator used to describe the linearized system dynamics. Clearly, for fluid mechanical systems, which in general present numerous degrees of freedom, the usage of a ROM is preferable (Kim & Bewley Reference Kim and Bewley2007).
For the realization of the ROM we use the ERA-POD (Juang & Pappa Reference Juang and Pappa1985). The ERA-POD is based on output signals resulting from impulse responses. It is necessary to have access to an number of impulse responses equal to the number of total inputs of the systems. The signals are written in a Hankel matrix whose dimensions depend on the total number of inputs and outputs and on the length of the saved time series, i.e. 
$N_{t}(N_{y}+N_{z})\times N_{t}(N_{u}+N_{d})$, 
$N_{t}$ being the number of time samples needed to have a good representation of the impulse response. In case the disturbance is free-stream turbulence the number of degrees of freedom used for the implementation of 
$\boldsymbol{d}(t)$ in the fully nonlinear N–S solver is of the order of hundreds. Then, since the Hankel matrix is decomposed by the singular value decomposition, it is clear that collecting such a high number of impulse responses results in a heavy computational problem. Besides, in a practical application it is not possible to collect the impulse responses from the free-stream turbulence disturbance.
Therefore, in order to reduce the computational power required and to have a method that can be applied in experiments as well, a different approach is proposed. A new set 
$N_{y}\times 1$ of outputs 
$\boldsymbol{y}_{d}(t)$, which measure the same quantity as 
$\boldsymbol{y}(t)$ and 
$\boldsymbol{z}(t)$, is introduced upstream of 
$\boldsymbol{y}(t)$, and the outputs generated by a nonlinear N–S simulation with free-stream turbulence and without control action are stored. An impulse response coincides with a TF by definition, so the following TFs can be computed as in (5.1),
$$\begin{eqnarray}\displaystyle & \displaystyle {\hat{G}}^{y_{d}y}=\frac{{\hat{S}}_{y_{d}y}}{{\hat{S}}_{y_{d}y_{d}}},\quad {\hat{G}}^{y_{d}z}=\frac{{\hat{S}}_{y_{d}z}}{{\hat{S}}_{y_{d}y_{d}}}, & \displaystyle\end{eqnarray}$$ and their inverse Fourier transforms can be used as a set of impulse responses to mimic the presence of free-stream turbulence upstream of every control device. These estimated impulse responses are used in the ERA-POD to model the impulse responses coming from 
$\unicode[STIX]{x1D648}_{d}\boldsymbol{d}(t)$ in (2.4), which represents the disturbance in the system. The number of impulse responses from the actuators, which are characterized by 
$\unicode[STIX]{x1D63D}\boldsymbol{u}(t)$, is 
$N_{u}$. Nevertheless, only one of these impulse responses is collected because the other ones can be computed by exploiting the homogeneity of 
$\boldsymbol{q}_{B}$ and the periodicity of the flow along the spanwise direction.
Once the whole set of impulse responses is available it is possible to build the mentioned Hankel matrix, apply the ERA-POD and retrieve the ROM needed for the design of the LQG. To the best of the authors’ knowledge, this is the first time this approach has been used in fluid mechanics.
(a–c) Comparison between empirical and improved TF; (a) empirical TF, 
$|{\hat{G}}_{emp}^{yz}|$; (b) improved TF, 
$|{\hat{G}}_{impr}^{yz}|$; (c) 
$|{\hat{G}}_{emp}^{yz}-{\hat{G}}_{impr}^{yz}|$. (d–f) Comparison between the true 
$z(t)_{k}$ output from DNS data at 
$x_{1}=400$ and the estimated output 
$\tilde{z}(t)_{k}$ (available output 
$y_{avail}(t)_{k}$ at 
$x_{1}=250$); (d) single output 
$k=9$; (e) empirical TF error (%), 
$||\hat{\tilde{z}}|^{2}-|\hat{z}|^{2}|/\max (S_{zz})$; (f) improved TF error (%), 
$||\hat{\tilde{z}}|^{2}-|\hat{z}|^{2}|/\max (S_{zz})$.
5.4 Identified reduced-order model
Given the position of the sensors, data-driven TFs can be computed by exploiting the methods outlined in § 5.1 and § 5.3 and improved as described in § 5.2. Ensemble averaging is performed on time series with a sampling frequency of 
$1/0.3$, over a sampling time of 
$T=25\,000$, with 16 000 samples per segment, an overlap of 
$80\,\%$, a total number of 22 segments and by means of a triangular windowing function. The improvement of the empirical TF, with available outputs 
$y(t)_{k}$ at 
$x_{1}=250$ and estimated outputs 
$\tilde{z}(t)_{k}$ at 
$x_{1}=400$, is summarized in figure 2 and table 1. Figure 2 shows a comparison between the TFs 
${\hat{G}}^{yz}$ estimated with the two methods and between the resulting estimated signals. It appears that the improved TF increases the accuracy in estimating the frequencies 
$(\unicode[STIX]{x1D714},\unicode[STIX]{x1D6FD}_{k})>(0.02,0.9)$ present in the original signal, as can be seen in figure 2(a–c). The reason is likely spectral leakage together with the lower amplitude that the higher frequencies have with respect to the lower frequencies. The error is reduced in the improved TF because it is computed in the physical domain, as in (5.2), where it is possible to isolate the erroneous frequencies bypassing the issue of relative amplitude and spectral leakage. Figure 2(d–f) shows that the signal estimated by the improved TF is closer to the measured signal, and highlights the frequencies where the error is reduced. It is noticeable that the error per frequency at 
$(\unicode[STIX]{x1D714},\unicode[STIX]{x1D6FD}_{k})<(0.02,0.9)$ is not reduced after the improvement. This suggests that the error at these frequencies belongs to the null space of the measurement operator (i.e. 
$\hat{p}$ in (3.5)), so it cannot be further reduced as discussed in § 3.2. Errors in amplitude estimation are also related to the windowing procedure of the time signal in the ensemble averaging. Even though the window is chosen to minimize such errors, its usage inevitably alters the computed amplitudes. Table 1 shows the normalized correlation value at zero delay, the mean square (
$\text{MS}$) and the variance (
$\text{VAR}$) for the empirical TF and for the improved TF.
The TF that estimates the output 
$\tilde{z}(t)_{k}$ given the available output 
$y(t)_{k}$, i.e. 
${\hat{G}}^{yz}$, is used in the IFFC control technique. The identified TFs used to mimic the effect of free-stream turbulence, characterized by 
$\unicode[STIX]{x1D648}_{d}\boldsymbol{d}(t)$ in the ROM, assume as available output a set of sensors 
$y_{d}(t)_{k}$ at 
$x_{1}=175$ and estimate the outputs at 
$x_{1}=250,400$, i.e. 
$\boldsymbol{y}(t)$ and 
$\boldsymbol{z}(t)$. These TFs correspond to (5.1).
The ROM resulting from the ERA-POD consists of 
$N=387$ degrees of freedom, which is considerably less than the degrees of freedom of the full system. The solution of the algebraic Riccati equations, which is the most computationally demanding step in the control design, with 
$N=387$ can be computed within the order of minutes nowadays (on a laptop). The value 
$N=387$ is found by imposing the error between the impulse response from the ROM and the original impulse response to be small enough. Since the ERA-POD performs the singular value decomposition of a Hankel matrix and the singular values are ordered such that 
$\unicode[STIX]{x1D70E}_{i}>\unicode[STIX]{x1D70E}_{i+1}$, the ratio 
$\unicode[STIX]{x1D70E}_{N}/\unicode[STIX]{x1D70E}_{max}\leqslant 5\times 10^{-4}$ is used to determine 
$N$. Moreover, given the equivalency between ERA-POD and approximate balanced POD truncation (Ma et al. Reference Ma, Ahuja and Rowley2011), the infinity norm of the error between the exact linear system and the ROM has the upper bound 
$2\sum _{i=1}^{\infty }\unicode[STIX]{x1D70E}_{N+i}\approx 0.56\times 10^{-2}\sum _{i=1}^{\infty }\unicode[STIX]{x1D70E}_{i}$ for the chosen 
$N$. Figure 3 shows the singular values of the Hankel matrix used for the ROM. Figure 4 compares the impulse response from 
$\boldsymbol{d}(t)$ in the ROM resulting from the ERA-POD against the estimated TF used as the original impulse response in the ERA-POD. The TFs are centred at zero and present a peak which is related to the group velocity of the structures. There clearly is good agreement between the ROM and the original data.
Singular values of the SVD of the Hankel matrix used for the ERA-POD. Dots, singular values. Dashed line, last singular value used for the ROM, 
$N=387$.
Correlation coefficient at zero delay (corresponding to its maximum value), mean square 
$\text{MS}[\boldsymbol{z}(t)]$ (DNS) or 
$\text{MS}[\tilde{\boldsymbol{z}}(t)]$ (estimation), and variance 
$\text{VAR}[\boldsymbol{z}(t)]$ (DNS) or 
$\text{VAR}[\tilde{\boldsymbol{z}}(t)]$ (estimation) for the cases shown in figure 2.
Original identified impulse responses used for the ERA-POD (solid black lines) versus ROM impulse response (dashed red lines). Impulse response from 
$\boldsymbol{d}(t)$ to 
$\boldsymbol{y}(t)$ and 
$\boldsymbol{z}(t)$. The original identified impulse responses are built by the improved frequency response technique. (a) Central output. (b) Complete set of TF.
6 Control performance: transition delay
Here, the results from the nonlinear simulations are presented. The controllers are designed on the ROM and then applied to the nonlinear N–S equations in a reduce-then-control approach (Anderson & Liu Reference Anderson and Liu1989). This approach can lead to inconsistencies because the disregarded modes of the full-order linearized system may become important in a closed loop, or the non-orthogonality of the linearized operator may lead to non-negligible nonlinearities in the closed-loop dynamics. Nevertheless, a Blasius boundary layer flow is convectively unstable and globally stable (Åkervik et al. Reference Åkervik, Ehrenstein, Gallaire and Henningson2008), so the full-order linearized system has all eigenvalues with negative real part, and the feed-forward configuration of the present plant grants the stability of the system as long as the controller is stable (Belson et al. Reference Belson, Semeraro, Rowley and Henningson2013; Semeraro et al. Reference Semeraro, Pralits, Rowley and Henningson2013b), which is assured by the control laws adopted. This implies that in the present flow case the reduce-then-control approach can lead to inconsistencies only if the closed-loop dynamics triggers non-negligible nonlinearities, which is shown not to be the case in the following sections. In fact, disturbance attenuation and transition delay are achieved. It should be noted that these non-negligible nonlinearities are strictly related to two factors: the spatial support of the actuator and the amplitude of actuation. An actuator can damp a disturbance in the linearized system, but requires an amplitude of actuation high enough to trigger nonlinearities in the original system. In this case the actuator is not efficient. This lack of efficiency lies in the spatial support of the actuator. In fact, it is the shape of the actuator that defines the controllability of the system (Lewis & Syrmos Reference Lewis and Syrmos1995; Bagheri et al. Reference Bagheri, Brandt and Henningson2009), which directly affects the gains of the controller that give the amplitude of actuation. Thus, the shape of the actuator is key to obtaining the desired performance. Here, the chosen actuator is efficient enough to require low amplitudes of the actuation signal, but care should be taken in the choice of the actuator (see Sasaki et al. (Reference Sasaki, Morra, Cavalieri, Hanifi and Henningson2019) for a discussion on the subject).
6.1 Transition delay
The 
$G_{m}^{yu}(t)$ used in the nonlinear N–S simulations are those resulting in the best performance in the control design. The weight matrices of the 
$J$ functionals of equations (3.4) and (3.8) are summarized in table 2 for IFFC and LQG, respectively. The choice of weights for control design is made through input–output simulations based on time signals stored from uncontrolled nonlinear N–S simulations. The input–output simulations consists of a linear superposition of signals, which makes them computationally inexpensive, such that their usage for control design becomes convenient to determine appropriate weights (details can be found in appendix C).
Design weights of the 
$G_{m}^{yu}(t)$ used in the fully nonlinear N–S simulations.
In order to assess the performance of the controller, the following quantity is introduced
$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{E}}=\frac{J_{controlled}^{M}}{J_{uncontrolled}^{M}},\quad J^{M}=\frac{1}{T}\int _{0}^{T}\boldsymbol{z}(t)^{T}\boldsymbol{z}(t)\,\text{d}t, & \displaystyle\end{eqnarray}$$ which corresponds to the average behaviour of the output to minimize. 
$T$ is the total time of the simulation.
Performance, 
${\mathcal{E}}$, for the two control methods. Input–output versus nonlinear N–S simulations.
Table 3 shows the performance of each control technique resulting from the input–output and nonlinear N–S simulations. The comparison of the control techniques in the input–output simulation is consistent with the results of the fully nonlinear N–S: LQG outperforms IFFC. Moreover, both cases present a smaller reduction of the objective function in the results from the N–S simulations than the input–output ones, which may be attributed to nonlinearities. In figure 5 the kernels from the IFFC and the LQG methods, respectively, are shown. These correspond to 
$G_{m}^{yu}(t)$ as in § 3. It appears that the LQG weights a lot more the recent history of the signal than the IFFC does, which may be seen as the main reason for outperforming the IFFC result, as further discussed next in § 6.2. In figure 6, the spanwise root-mean-square (r.m.s.) values of 
$G_{m}^{yu}(t)$
$$\begin{eqnarray}\displaystyle & \displaystyle \text{rms}[u]_{x_{3}}=\left(\frac{1}{N_{u}}\mathop{\sum }_{k=1}^{N_{u}}u(t)_{k}^{2}\right)^{1/2}=\left(\frac{1}{N_{u}}\boldsymbol{u}(t)^{T}\boldsymbol{u}(t)\right)^{1/2}, & \displaystyle\end{eqnarray}$$ are plotted. We observe that: (i) the magnitude of the averages and the fluctuation around the average values are higher for the signal generated by 
$G_{m}^{yu}(t)$ from the LQG, and (ii) the actuation signal from the IFFC is smoother. The first difference can explain why the performance of the LQG drops more than that of the IFFC when moving from the input–output to the N–S simulations. In fact, the actuation signal multiplies a fixed spatial support, and an increase in the magnitude of the actuation signal corresponds to a more intense forcing that leads to stronger nonlinearities. The second difference comes from the shape of the 
$G_{m}^{yu}(t)$ in figure 5. The 
$G_{m}^{yu}(t)$ from the IFFC is less localized around 
$t=0$ than the one from the LQG, which means that the latter only captures the low-frequency dynamics of 
$y(t)_{k}$ signals. This can also be seen by comparing figures 6 and 7.
Control TF 
$G_{k}^{yu}(t)$. (a) IFFC technique. (b) LQG technique. Notice the different scales for the gains between the two panels.
Spanwise r.m.s. values of the actuation signals 
$u(t)_{k}$ at 
$x_{1}=325$. Solid line, LQG. Dotted line, IFFC.
Spanwise r.m.s. values of the output 
$y(t)_{k}$ at 
$x_{1}=250$.
(a) 
$q_{1,rms}^{\prime }/{q^{\prime }}_{1,rms}^{uncontrolled}$ averaged along the spanwise direction at 
$Re_{x}=(1.51,1.96,2.40,2.86)\times 10^{5}$; red dash-dotted line, uncontrolled case; black solid line, LQG. (b) 
$q_{1,rms}$ averaged along the spanwise direction at 
$Re_{x}=1.5\times 10^{5}$; red dash-dotted line, uncontrolled case; black dotted line, IFFC; black solid line, LQG.
The effect of the actuation signal on the spanwise r.m.s. values of the streamwise velocity component
$$\begin{eqnarray}\displaystyle q_{1,rms}^{\prime }=\sqrt{\langle (q_{1}^{\prime }-\langle q_{1}^{\prime }\rangle _{t})^{2}\rangle _{t,x_{3}}}, & & \displaystyle\end{eqnarray}$$ with 
$\langle \bullet \rangle$ representing the sample average, at 
$Re_{x}=(1.51,1.96,2.40,2.86)\times 10^{5}$, is shown in figure 8. It is clear that LQG outperforms slightly IFFC also in terms of reduction of the disturbance amplitude throughout the boundary layer. The disturbance energy drops by a factor of 
${\approx}40\,\%$ after the control action. As shown in figure 9, despite the growth of disturbance amplitude beyond 
$x_{1}=400$ (
$Re_{x}\approx 1.5\times 10^{5}$) where the output 
$z(t)_{k}$ for the objective function is measured, the smaller amplitude of the streamwise component of the perturbation and the downstream shift of the curves imply dampening of the disturbance, and thus transition delay.
Maximum along the wall-normal direction of the 
$q_{1,rms}$ averaged along the spanwise direction. Red dash-dotted line, uncontrolled case. Black dotted line, IFFC. Black solid line, LQG.
A different measure of transition delay is the skin friction coefficient, which explicitly appears in the computation of the drag and is the measure of interest for many applications. A measure of the skin friction coefficient based on the r.m.s. values of the streamwise velocity component
$$\begin{eqnarray}\displaystyle & \displaystyle C_{f}(x_{1})=\left[\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{2}}\sqrt{\langle {q^{\prime }}_{1}^{2}\rangle _{t,x_{3}}}\right]_{x_{2}=0} & \displaystyle\end{eqnarray}$$ is shown in figure 10. There, the threshold curves that represent the skin friction of a laminar and a fully turbulent flat-plate boundary layer are also presented. It is clearly the case again that LQG performs slightly better than IFFC and that in the best case the transition delay is around 
$\unicode[STIX]{x0394}Re_{x}=1.5\times 10^{5}$, which is at least as good as most of the current results in the literature where more idealized cases are studied, as in Monokrousos et al. (Reference Monokrousos, Brandt, Schlatter and Henningson2008). There, a transition delay around 
$\unicode[STIX]{x0394}Re_{x}=1.2\times 10^{5}$ is achieved in the best possible scenario for a case with turbulence intensity 
$Tu=3.0\,\%$ and integral length scale 
$L=5.0\unicode[STIX]{x1D6FF}_{0}^{\ast }$. Nevertheless, in Monokrousos et al. (Reference Monokrousos, Brandt, Schlatter and Henningson2008), actuation is performed by control of each point on a band of the flat plate, while measurements consists of wall shear stress in both streamwise and spanwise directions and wall pressure fluctuations over a band of the flat plate located downstream of the actuation.
Average skin friction coefficient 
$C_{f}$. Red dash-dotted line, uncontrolled case. Black dotted line, IFFC control method. Black solid line, LQG control method.
At a fixed 
$Re_{x}$ the mean value of the 
$C_{f}$ of a fully turbulent and a laminar boundary layer can be taken as reference to indicate a significant appearance of turbulence spots. This mean value is reached at 
$Re_{x}\approx 5\times 10^{5}$ for the uncontrolled case and at 
$Re_{x}\approx 6.5\times 10^{5}$ for the controlled case with LQG, as shown in figure 10. These values represent only an indication of the significant appearance of turbulence spots. The instantaneous behaviour of the skin friction coefficient over the flat plate
$$\begin{eqnarray}\displaystyle & \displaystyle c_{f}(x_{1},x_{3},t)=\left[\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{2}}q_{1}^{\prime }(\boldsymbol{x},t)\right]_{x_{2}=0} & \displaystyle\end{eqnarray}$$ is shown in figure 11. The uncontrolled flow shows turbulent spots upstream of 
$Re_{x}=5\times 10^{5}$ together with wiggles typically present when secondary instability occurs, as expected in laminar-to-turbulent transition due to the breakdown of streaky structures. Only chaotic structures appear for 
$Re_{x}>5\times 10^{5}$, consistently with the trend of 
$C_{f}$ in figure 10. The controlled flow shows the same behaviour of the uncontrolled one but further downstream, with chaotic structures appearing only in an upstream neighbourhood of 
$Re_{x}=6.5\times 10^{5}$. The controlled flow clearly appears smoother than the uncontrolled one for 
$Re_{x}<6\times 10^{5}$, without the appearance of secondary instabilities up to 
$Re_{x}\approx 5.2\times 10^{5}$. This is caused by the lower streak amplitude (less pronounced white and black areas in figure 11) achieved thanks to the application of the controller.
Instantaneous skin friction coefficient 
$c_{f}(Re_{x},x_{3},t=t^{\ast })$. (a) Uncontrolled case. (b) Controlled case with LQG. Black and white colours, 
$c_{f}<2.5\times 10^{-3}$. Empty contours, 
$c_{f}\geqslant 2.5\times 10^{-3}$. Same time 
$t=t^{\ast }$ for both simulations; the starting seed for the random free-stream turbulence generation is the same.
6.2 Role of control methods, sensors and actuators for control performance
Figure 12 shows the difference between the uncontrolled and the controlled fields, 
$q_{1,unc}^{\prime }-q_{1,ctr}^{\prime }$, at a time 
$t=t^{\ast }$ for the LQG and the IFFC cases. Such fields can be understood as the disturbances induced by the actuators in the simulations. Both controllers appear to trigger streaky structures in the boundary layer flow. Some differences between the LQG and the IFFC actions appear, but cannot be attributed to the stochastic nature of the disturbance because the starting seed for the random free-stream turbulence is the same for all simulations. The actuator is the same in both cases and the controllers handle only their amplitudes, so the differences between the actions of the LQG and the IFFC must reside in the time signals. These differences must also explain the reason behind the LQG outperforming the IFFC, and is now discussed from a physical point of view.
A limiting characteristic for control performance is the relative position of sensors and actuators, which was chosen to minimize prediction errors and exploit the linear behaviour of the flow field. The input-to-output and the output-to-output time delays are respectively 
$\unicode[STIX]{x1D70F}_{uz}=219.6$ and 
$\unicode[STIX]{x1D70F}_{yz}^{IFFC}=216$; these delays correspond to the peak of the transfer function, and approximate the average time for a structure induced by free-stream turbulence or the actuator, respectively, to arrive at the 
$\boldsymbol{z}(t)$ sensors. Therefore, the actuation is not sufficiently fast to cancel the streaks detected by the 
$\boldsymbol{y}(t)$ sensors, which is reflected by IFFC resulting in a non-causal 
$G^{yu}$, as shown in figure 13. The result suggests that the performance of the controller may improve if it were possible to increase the difference between the time delays, which can be achieved either by changing the relative position of the devices or by changing the type of sensors or actuators. The limitations in the control performance are therefore not caused by the control methods, but by the structure of the plant. However, it appears that LQG can slightly compensate for this causality issue without any modification to the plant.
Effect of the actuator. Instantaneous subtraction between an uncontrolled and a controlled field, 
$q_{1,unc}^{\prime }(\boldsymbol{x},t^{\ast })-q_{1,ctr}^{\prime }(\boldsymbol{x},t^{\ast })$, at a time 
$t=t^{\ast }$ for the LQG (a,c) and the IFFC (b,d). Streamwise component of the perturbation velocity. (a,b) Cross-plane at 
$x_{1}=400\unicode[STIX]{x1D6FF}_{0}^{\ast }$. (c,d) Wall-parallel plane at 
$x_{2}=2.0\unicode[STIX]{x1D6FF}_{0}^{\ast }$. Red and solid lines, positive values; blue and dashed lines, negative values. The starting seed for the random free-stream turbulence generation is the same.
Full non-causal control TF 
$G_{k}^{yu}(t)$ from the IFFC control method. (a) Central line; the dot represents the peak value. (b) Full TF.
There exists a specific set of weights for which LQG outperforms IFFC, but there also exists a different set of weights for which LQG results in the same solution given by the IFFC. This is possible thanks to the presence of the estimation problem in the LQG, which introduces two more degrees of freedom in the control design. It follows that the reason behind the better performance of the LQG is the possibility of optimizing the estimation inside the control method, which is not included in the IFFC. Moreover, it appears that in the LQG keeping fixed the weight ratio 
$R^{LQG}/Q^{LQG}$ associated with the best solution and increasing the value of the ratio 
$v_{n}/v_{d}$ leads to worse performance (appendix C). Thus, since the weights of the estimation problem define the dynamics of the estimator (appendix A), the best solution comes from the estimator with the fastest response. In the present ROM the estimator that corresponds to 
$\unicode[STIX]{x1D70F}_{yz}^{LQG}=\unicode[STIX]{x1D70F}_{yz}^{IFFC}=216$, as in the 
$G^{yz}$ used in IFFC, has a slower response than the one corresponding to the best solution found with LQG, where 
$\unicode[STIX]{x1D70F}_{yz}^{LQG}=213$. Thus, it appears that LQG achieves better performance for a case where the difference between the time delays, 
$\unicode[STIX]{x1D70F}_{uz}-\unicode[STIX]{x1D70F}_{yz}$, is higher than it is in IFFC, as suggested by the non-causal result found from the IFFC, and that better performance is achieved by an estimator with a fast response. The latter occurs because the weights of the estimation problem define the dynamics of its error: an estimator with a fast response has a fast decaying error. This implies that after the same 
$\unicode[STIX]{x0394}t$ the faster estimator is more accurate.
(a) Cross-correlation between the streamwise velocity fluctuation at the peak of the 
$u_{rms}$ profile and the shear of the streamwise fluctuation at the wall. (b) Contours are the fluctuations of the streamwise velocity component; red and solid lines, positive values; blue and dashed lines, negative values; dashed-dotted line, 
$\unicode[STIX]{x1D6FF}_{99}/\unicode[STIX]{x1D6FF}_{0}^{\ast }$, boundary layer.
The connection between the shorter time delay 
$\unicode[STIX]{x1D70F}_{yz}^{LQG}=213$ mentioned earlier and the presence of a fast responding estimator may be explained by analysing the capability of the wall-shear-stress measurements to capture the dynamics of the streaks. An alternative measure of the streaky perturbations is their maximum streamwise velocity. Thus, a new set of outputs at 
$(x_{1},x_{2})=(250,2.25)$ with the same spanwise positions of the other considered outputs is collected, as in figure 14(b). These outputs are compared to those that measure the shear on the wall at 
$(x_{1},x_{2})=(250,0)$, which are used to compute the input 
$u(t)_{k}$. Figure 14(a) shows the time–space correlation between the two sets of measurements. It appears that the output resulting from the measurements of the streamwise velocity perturbation and the output resulting from the measurement of the shear on the wall are highly correlated in the positive time half-plane. With the considered convention for the correlation, this implies that fluctuations at the higher wall-normal position precede those on the wall, and thus the instantaneous measure of the shear on the wall cannot predict the instantaneous or future maximum velocity fluctuation of the streak. In other terms, there is reverse causality between the measurements at 
$x_{2}=0$ and the measurements at 
$x_{2}=2.25$. This effect can be associated with the tilting of the streaky structures shown in figure 14(b): an advecting streak first passes at higher wall-normal positions (exemplified by the considered probe), and only later leaves a wall-shear-stress signature. However, the wall-shear-stress measurements are not completely unable to estimate the streaks’ velocity; they can effectively predict the velocity of the tilted structure at wall-normal positions closer to the wall. There the velocity is lower, the convection velocity of the streaks is reduced and thus the time delay 
$\unicode[STIX]{x1D70F}_{yz}$, which describes their travelling time, is larger. A more accurate estimation, which corresponds to an estimator with a faster response, can effectively predict the velocity further from the wall. There the velocity is higher, so closer to the real travelling speed of the streaks, and the resulting time delay 
$\unicode[STIX]{x1D70F}_{yz}$ is smaller. This explains why to the best LQG solution corresponds a fast estimator and its connection to a shorter time delay 
$\unicode[STIX]{x1D70F}_{yz}^{LQG}=213$.
Moreover, the fact that the LQG results in a 
$G^{yu}$ which puts a lot of weight onto the recent history of the output signal, around two orders of magnitude higher than the one from the IFFC (figure 5), is also explained by the presence of a fast estimator. In fact, the higher values of 
$v_{n}/v_{d}$, with 
$R_{LQG}/Q_{LQG}$ fixed to the value of the best solution, correspond to a shape of 
$G^{yu}$ which approaches the one of the best IFFC solution, so the difference between the best LQG and IFFC results must come from the presence of the fast estimator. The last statement is consistent with the present discussion, as it implies that the performance of the controller improves when it mainly makes use of the portion of the output that lies in a small neighbourhood of 
$t=T$, with 
$t\in (-\infty ,T]$ the history time and 
$T$ the running time. This neighbourhood contains the meaningful information because of the mentioned reverse causality between the wall-shear-stress measurements and the dynamics of the streaks, as shown by the non-positive time half-plane in figure 14(a).
In summary, it appears that using streamwise wall shear stress as a measure to predict the dynamics of the streaks introduces a time delay in the ROM because of the tilting of the structures. The time delay affects the actuation signal because the controller is designed on the ROM, so the streaky structures generated by the actuators are not exactly in phase with the incoming disturbance. Moreover, it can be concluded that the LQG can achieve better performance than the IFFC thanks to the weights in the design of the estimator. By choosing the appropriate value for these weights the estimation accuracy increases because the estimator slightly compensates for the reverse causality between the wall-shear-stress measurements and the dynamics of the streaks. Finally, it appears that the limitation caused by the structure of the plant, including the location of sensors and actuators and their shapes, is more critical than the choice of control technique, and thus is the key design challenge (as further discussed in the parallel work of Sasaki et al. (Reference Sasaki, Morra, Cavalieri, Hanifi and Henningson2019)).
7 Conclusions
The delay of bypass transition in a realistic scenario by means of active flow control, control theory and system identification is presented. Numerical simulations of the nonlinear transitional regime in a Blasius boundary layer are performed, where streaks are excited in the boundary layer by means of a high level of free-stream turbulence. A model-based method for the delay of bypass transition realizable in experiments is introduced. It makes use of a ROM representation of the system and is based on the signals from a finite number of localized sensors and actuators placed on the wall, which mimic real shear-stress sensors and ring plasma actuators. A technique for the characterization of disturbances with a large number of degrees of freedom for model-based approaches is presented, which allows us to obtain reasonably low-dimensional ROMs by isolating the dynamics of interest via system identification of the effects of the disturbance on the system. The method is reliable, easy to implement and based on measurement data, in a data-driven approach that would be realizable in experiments. The presented technique is applied to generate the ROM via ERA-POD for solving the flow control problem by means of LQG, which to the best of the authors’ knowledge has never been done in a flow control application.
The performance of the LQG is compared to that of the IFFC optimal control technique, which does not need the explicit characterization of the disturbance on the system, thus, simplifying the flow control problem. The LQG is seen to perform slightly better than the simpler IFFC method once appropriate weights in the cost function are selected. The performance of the control techniques are compared in linear input–output and nonlinear N–S simulations, showing that resorting to a linear ROM for control design is reasonable also in presence of the high-amplitude disturbances considered here. The effectiveness of the technique in delaying bypass transition is shown. Using LQG a transition delay of 
$\unicode[STIX]{x0394}Re_{x}\approx 1.4\times 10^{5}$ for a case with turbulent intensity 
$Tu=3.0\,\%$ and integral length scale 
$L=7.5\unicode[STIX]{x1D6FF}_{0}^{\ast }$ is achieved. This highlights the capability of the presented methods to achieve at least as large a delay of bypass transition as that obtained in more idealized cases found in the literature (Monokrousos et al. Reference Monokrousos, Brandt, Schlatter and Henningson2008).
Finally, the differences in the results obtained with IFFC and LQG are analysed and related to the structure of the plant, so the limitations caused by the relative positions of sensors and actuators and by the shape of the sensor are outlined. In particular, a reverse causality issue arising from using wall streamwise-shear-stress sensors to predict the dynamics of the streaks is shown. The way in which this causality issue limits the control performance is described and an explanation on the way in which the LQG can compensate for such issue is provided.
Acknowledgements
The authors would like to acknowledge the VINNOVA Projects PreLaFlowDes and SWE-DEMO and the Swedish-Brazilian Research and Innovation Centre CISB for funding. Moreover, part of this work was performed during an exchange programme at KTH, for which K.S. received a scholarship from Capes, project number 88881.132008/ 2016-01. K.S. work is also funded by a scholarship from FAPESP, grant number 2016/25187-4. A.V.G.C. was supported by a CNPq grant 310523/2017-6. The simulations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at NSC, HPC2N and PDC.
Appendix A. Linear quadratic Gaussian regulator
The LQG technique is designed to solve the control problem on a dynamical system subject to stochastic white noise disturbances. Here, the dynamical system is a ROM and reads
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \dot{\boldsymbol{q}}(t)=\unicode[STIX]{x1D63C}\boldsymbol{q}(t)+\unicode[STIX]{x1D63D}\boldsymbol{u}(t)+\unicode[STIX]{x1D648}_{d}\boldsymbol{d}(t),\\ \displaystyle \boldsymbol{y}(t)=\unicode[STIX]{x1D63E}_{y}\boldsymbol{q}(t)+\boldsymbol{n}(t),\\ \displaystyle \boldsymbol{z}(t)=\unicode[STIX]{x1D63E}_{z}\boldsymbol{q}(t).\end{array}\right\} & & \displaystyle\end{eqnarray}$$Since the LQG does not assume the full state to be known, an estimation of the original dynamical system based on the known outputs is introduced,
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \dot{\tilde{\boldsymbol{q}}}(t)=\unicode[STIX]{x1D63C}\tilde{\boldsymbol{q}}(t)+\unicode[STIX]{x1D63D}\boldsymbol{u}(t)-\unicode[STIX]{x1D647}(\,\boldsymbol{y}(t)-\tilde{\boldsymbol{y}}(t)),\\ \displaystyle \tilde{\boldsymbol{y}}(t)=\unicode[STIX]{x1D63E}_{y}\tilde{\boldsymbol{q}}(t),\\ \displaystyle \tilde{\boldsymbol{z}}(t)=\unicode[STIX]{x1D63E}_{z}\tilde{\boldsymbol{q}}(t),\end{array}\right\} & & \displaystyle\end{eqnarray}$$ where 
$\tilde{\boldsymbol{q}}(t)$ and 
$\tilde{\boldsymbol{y}}(t)$ are estimates of 
$\boldsymbol{q}(t)$ and 
$\boldsymbol{y}(t)$, and 
$\unicode[STIX]{x1D647}$ is an 
$N\times N_{y}$ matrix to be designed. The estimated system accounts for the stochastic disturbances through the available outputs 
$\boldsymbol{y}(t)$. Subtracting (A 2) from (A 1) and substituting 
$\boldsymbol{y}(t)=\unicode[STIX]{x1D63E}_{y}\boldsymbol{q}(t)$ and 
$\tilde{\boldsymbol{y}}(t)=\unicode[STIX]{x1D63E}_{y}\tilde{\boldsymbol{q}}(t)$ gives
$$\begin{eqnarray}\displaystyle \dot{\boldsymbol{e}}(t)=(\unicode[STIX]{x1D63C}+\unicode[STIX]{x1D647}\unicode[STIX]{x1D63E}_{y})\boldsymbol{e}(t)+\unicode[STIX]{x1D648}_{d}\boldsymbol{d}(t)+\unicode[STIX]{x1D647}\boldsymbol{n}(t), & & \displaystyle\end{eqnarray}$$ with 
$\boldsymbol{e}(t)=\boldsymbol{q}(t)-\tilde{\boldsymbol{q}}(t)$ the estimation error. Equation (A 3) shows that the error dynamics is based on the matrix 
$\unicode[STIX]{x1D647}$ and is driven by the stochastic disturbances. Thus, the matrix 
$\unicode[STIX]{x1D647}$ should stabilize the error dynamics and dampen the amplitude of the stochastic disturbance 
$\boldsymbol{n}(t)$.
Since 
$\boldsymbol{y}(t)$ is an available measure, the estimated system is deterministic. Its solution is used to compute the actuation input 
$\boldsymbol{u}(t)=\unicode[STIX]{x1D646}\tilde{\boldsymbol{q}}(t)$, with 
$\unicode[STIX]{x1D646}$ a matrix to be designed to solve the control problem. Substituting 
$\boldsymbol{u}(t)=\unicode[STIX]{x1D646}\tilde{\boldsymbol{q}}(t)$ in (A 2) gives
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \dot{\tilde{\boldsymbol{q}}}(t)=(\unicode[STIX]{x1D63C}+\unicode[STIX]{x1D63D}\unicode[STIX]{x1D646})\tilde{\boldsymbol{q}}(t)-\unicode[STIX]{x1D647}(\,\boldsymbol{y}(t)-\tilde{\boldsymbol{y}}(t)),\\ \displaystyle \tilde{\boldsymbol{y}}(t)=\unicode[STIX]{x1D63E}_{y}\tilde{\boldsymbol{q}}(t),\\ \displaystyle \tilde{\boldsymbol{z}}(t)=\unicode[STIX]{x1D63E}_{z}\tilde{\boldsymbol{q}}(t).\end{array}\right\} & & \displaystyle\end{eqnarray}$$ The LQG technique consists of computing 
$\unicode[STIX]{x1D646}$ and 
$\unicode[STIX]{x1D647}$ to solve the control and the estimation problem, respectively. These two problems are usually coupled in optimal control, but in the LQG technique they are not and they result in the minimization of two different 
${\mathcal{H}}_{2}$-norms (Skogestad & Postlethwaite Reference Skogestad and Postlethwaite2005). The matrix 
$\unicode[STIX]{x1D646}$ results from the linear quadratic regulator problem. It minimizes the objective function
$$\begin{eqnarray}\displaystyle & \displaystyle J=\lim _{T\rightarrow \infty }\frac{1}{T}\int _{0}^{T}(\boldsymbol{z}(t)^{T}\unicode[STIX]{x1D64C}\boldsymbol{z}(t)+\boldsymbol{u}(t)^{T}\unicode[STIX]{x1D64D}\boldsymbol{u}(t))\,\text{d}t, & \displaystyle\end{eqnarray}$$and results in solving the following algebraic Riccati equation
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D63C}^{\text{T}}\unicode[STIX]{x1D64B}_{u}+\unicode[STIX]{x1D64B}_{u}\unicode[STIX]{x1D63C}-\unicode[STIX]{x1D64B}_{u}\unicode[STIX]{x1D63D}\unicode[STIX]{x1D64D}^{-1}\unicode[STIX]{x1D63D}^{\text{T}}\unicode[STIX]{x1D64B}_{u}+\unicode[STIX]{x1D63E}_{z}^{\text{T}}\unicode[STIX]{x1D64C}\unicode[STIX]{x1D63E}_{z}=0, & \displaystyle\end{eqnarray}$$ where 
$\unicode[STIX]{x1D64B}_{u}$ is a positive semi-definite 
$N\times N$ matrix which is the unknown of the equation. The relationship between 
$\unicode[STIX]{x1D646}$ and 
$\unicode[STIX]{x1D64B}_{u}$ reads
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D646}=-\unicode[STIX]{x1D64D}^{-1}\unicode[STIX]{x1D63D}^{\text{T}}\unicode[STIX]{x1D64B}_{u}. & \displaystyle\end{eqnarray}$$ The matrix 
$\unicode[STIX]{x1D647}$ results from the Kalman filter. It minimizes the expected value of the covariance matrix of the error at steady state
$$\begin{eqnarray}\displaystyle & \displaystyle J=\lim _{t\rightarrow \infty }\text{Tr}(\unicode[STIX]{x1D64B}_{e}(t)),\quad \unicode[STIX]{x1D64B}_{e}=\mathbb{E}[\boldsymbol{e}(t)\boldsymbol{e}(t)^{\text{T}}], & \displaystyle\end{eqnarray}$$ with 
$\text{Tr}(\bullet )$ the trace operator, and results in solving the following algebraic Riccati equation
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D64B}_{e}\unicode[STIX]{x1D63C}^{\text{T}}+\unicode[STIX]{x1D63C}\unicode[STIX]{x1D64B}_{e}-\unicode[STIX]{x1D64B}_{e}\unicode[STIX]{x1D63E}_{y}^{\text{T}}\unicode[STIX]{x1D651}_{n}^{-1}\unicode[STIX]{x1D63E}_{y}\unicode[STIX]{x1D64B}_{e}+\unicode[STIX]{x1D648}_{d}\unicode[STIX]{x1D651}_{d}\unicode[STIX]{x1D648}_{d}^{\text{T}}=0, & & \displaystyle\end{eqnarray}$$ where 
$\unicode[STIX]{x1D64B}_{e}$ is a positive semi-definite 
$N\times N$ matrix and is the unknown of the equation, and 
$\unicode[STIX]{x1D651}_{d}$ and 
$\unicode[STIX]{x1D651}_{n}$ are the covariance matrices of 
$\boldsymbol{d}(t)$ and 
$\boldsymbol{n}(t)$, respectively. The relationship between 
$\unicode[STIX]{x1D647}$ and 
$\unicode[STIX]{x1D64B}_{e}$ reads
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D647}=-\unicode[STIX]{x1D64B}_{e}\unicode[STIX]{x1D63E}_{y}^{\text{T}}\unicode[STIX]{x1D651}_{n}^{-1}. & & \displaystyle\end{eqnarray}$$ Once both 
$\unicode[STIX]{x1D646}$ and 
$\unicode[STIX]{x1D647}$ are computed the state-space system based on 
$\tilde{\boldsymbol{q}}$ gives the input signal based on the history of the available output 
$\boldsymbol{y}(t)$,
$$\begin{eqnarray}\displaystyle \boldsymbol{u}(t)=-\int _{0}^{t}\unicode[STIX]{x1D646}\text{e}^{(\unicode[STIX]{x1D63C}+\unicode[STIX]{x1D63D}\unicode[STIX]{x1D646}+\unicode[STIX]{x1D647}\unicode[STIX]{x1D63E}_{y})(t-\unicode[STIX]{x1D70F})}\unicode[STIX]{x1D647}\boldsymbol{y}(\unicode[STIX]{x1D70F})\,\text{d}\unicode[STIX]{x1D70F}. & & \displaystyle\end{eqnarray}$$Appendix B. Prediction for different streamwise positions
Performance of the estimation of downstream outputs from available upstream outputs: 
$x_{1}^{up}$ represents the available output position; 
$x_{1}^{down}$ represents the position of the estimated output. The outputs are the shear on the wall.
The streamwise position of sensors and actuators on the flat plate was chosen by also taking into account the estimation error. The behaviour of the estimation error as a function of the relative streamwise position between the available output and the output to estimate and as a function of the absolute position of the available output is inspected. Given the set of 
$N_{y}=36$ estimated outputs 
$\tilde{\boldsymbol{y}}_{est}^{\star }(t)$ and the set of true outputs 
$\boldsymbol{y}_{DNS}^{\star }(t)$, the error is defined as 
$\unicode[STIX]{x0394}\boldsymbol{y}^{\star }(t)=\boldsymbol{y}_{est}^{\star }(t)-\boldsymbol{y}_{DNS}^{\star }(t)$. The available outputs are placed at a streamwise position 
$x_{1}=x_{1}^{up}$ and the estimated outputs 
$\boldsymbol{y}_{est}^{\star }(t)$ are placed at a streamwise position 
$x_{1}=x_{1}^{down}$. The true outputs 
$\boldsymbol{y}_{DNS}^{\star }(t)$ are at the same place as the available outputs. It also holds that 
$x_{1}^{down}>x_{1}^{up}$, so the position of the outputs to be estimated never coincides with that of the available outputs. The available outputs are used to predict the downstream outputs in the future. Figure 15 shows the 
$\text{MS}[\unicode[STIX]{x0394}\boldsymbol{y}^{\star }(t)]$. Estimation is performed by means of empirical TFs because of their low computational cost. In figure 15 it is evident that the current positioning of sensors and actuators, 
$250\leqslant x_{1}\leqslant 400$, is adequate. For 
$x_{1}^{up}=150,250$, the 
$\text{MS}[\unicode[STIX]{x0394}\boldsymbol{y}^{\star }(t)]$ initially decreases due to decay of free-stream turbulence intensity along the streamwise direction. For 
$x_{1}^{down}>350$, the 
$\text{MS}[\unicode[STIX]{x0394}\boldsymbol{y}^{\star }(t)]$ increases in all cases because from that position the nonlinear interactions become important. The value of 
$\text{MS}[\unicode[STIX]{x0394}\boldsymbol{y}^{\star }(t)]$ grows faster with increasing 
$x_{1}^{down}$ as the flow nonlinearity increases.
In order to further confirm the fact that the chosen positioning of sensors and actuators is adequate, the coherence coefficient between the measurements at 
$x_{1}=250$ and 
$x_{1}=400$ is calculated. The coherence 
$\unicode[STIX]{x1D6FE}_{yz}$ is defined as
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FE}_{yz}^{2}=\frac{|{\hat{S}}_{yz}|^{2}}{{\hat{S}}_{yy}{\hat{S}}_{zz}}. & \displaystyle\end{eqnarray}$$ The definition holds for any pair of outputs. Its value varies between zero and one and indicates that the TF may be nonlinear, noise enters the measurements or 
$\boldsymbol{z}$ is the result of other inputs as well.
It is desirable to have the highest values of coherence in 
$(\unicode[STIX]{x1D714},\unicode[STIX]{x1D6FD}_{k})$ regions where the signals are most energetic. This may be evaluated by computing the power spectral density (PSD). The coherence coefficient between signals at 
$x_{1}=250$ and 400 and its normalized PSD is presented in figure 16. The results indicate an almost linear relation between signals at these two streamwise positions for part of the 
$(\unicode[STIX]{x1D714},\unicode[STIX]{x1D6FD}_{k})$ space that is of interest. It is noticeable that most of the energy is strongly localized at very low frequencies, 
$\unicode[STIX]{x1D714}\approx 0$, with a spanwise wavenumber 
$\unicode[STIX]{x1D6FD}\approx 0.4$, which corresponds to streaky motions with slow streamwise variation. This is the reason for the observed good accuracy of the estimation and justifies the choices of placements for sensors and actuators.
(a) Coherence 
$\unicode[STIX]{x1D6FE}_{yz}$ between the outputs 
$y(t)_{k}$, at 
$x_{1}=250$ and 
$z(t)_{k}$, at 
$x_{1}=400$. (b) Normalized power spectral density of the output 
$z(t)_{k}$: 
${\hat{S}}_{zz}/\max ({\hat{S}}_{zz})$.
Appendix C. Choice of weights
The LQG has two objective functions, one for the control problem and one for the estimation problem, whereas the IFFC has only one objective function, the one for the control problem. The objective function for the control problem requires the definition of the weight matrix on the output 
$\boldsymbol{z}(t)$ (or 
$\hat{z}$), 
$\unicode[STIX]{x1D64C}$ (or 
$\hat{Q}$), and the penalization matrix on the input 
$\boldsymbol{u}(t)$ (or 
$\hat{u}$), 
$\unicode[STIX]{x1D64D}$ (or 
$\hat{R}$), while the estimation problem requires the matrices that describe the covariance of the stochastic disturbance 
$\boldsymbol{d}$, 
$\unicode[STIX]{x1D651}_{d}$, and the noise 
$\boldsymbol{n}(t)$, 
$\unicode[STIX]{x1D651}_{n}$. The control problem deals with finding the function that, given an output, provides an input to minimize an objective function, while the estimation problem deals with finding the function that allows us to minimize the error in the estimation. The weights introduce more degrees of freedom in the design problem, and are usually left as free parameters. In fact, there is not a universally acclaimed method to compute those weights and close the control design problem, such that they are usually chosen iteratively (Skogestad & Postlethwaite Reference Skogestad and Postlethwaite2005). Here, a brute force method is applied: a grid of arbitrarily chosen weights is used and a set of 
$G_{m}^{yu}(t)$ is computed by means of the two control techniques as in § 3. The computed 
$G_{m}^{yu}(t)$ are tested in input–output simulations based on the linear superposition of the input and output time series only. The input–output simulations make use of the second in (2.7) to compute the effect of 
$G_{m}^{yu}(t)$ on the reference output to annihilate, 
$\boldsymbol{z}(t)$, and on (3.2) for the relationship between 
$\boldsymbol{u}(t)$ and 
$\boldsymbol{y}(t)$. The input–output simulation consists of computing for each time step
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle u(t)_{k}=\int _{0}^{t}\mathop{\sum }_{m=1}^{N_{y}}G_{m}^{yu}(t-\unicode[STIX]{x1D70F})y(\unicode[STIX]{x1D70F})_{m+k-1}\,\text{d}\unicode[STIX]{x1D70F},\\ \displaystyle z(t)_{k}=\int _{0}^{t}\mathop{\sum }_{m=1}^{N_{u}}G_{m}^{uz}(t-\unicode[STIX]{x1D70F})u(\unicode[STIX]{x1D70F})_{m+k-1}\,\text{d}\unicode[STIX]{x1D70F}+z^{d}(t)_{k},\end{array}\right\} & & \displaystyle\end{eqnarray}$$ with 
$y(t)_{k}$ and 
$z^{d}(t)_{k}$, at 
$x_{1}=250,400$, respectively, with the time series of outputs saved from the nonlinear uncontrolled N–S simulations. The method is thus a simpler simulation of the control effect considering only a linear superposition of the open-loop output 
$z^{d}(t)$ and what would result from control action (via the transfer function 
$G_{m}^{uz}$). The 
$k$ index was dropped in 
$G_{km}^{uz}(t)$ from (3.2) because the actuators have all the same spatial support and the linearized system dynamics is instantaneously homogeneous along the spanwise direction 
$x_{3}$. Here, 
$G_{m}^{uz}(t)$ is found from an impulse-response simulation of the linearized N–S equations.
This method avoids the use of computationally demanding N–S simulations and is reliable in identifying the best 
$G_{m}^{yu}(t)$ and the associated weights. It also proves to be consistent with the results of the nonlinear N–S simulations. The time required to perform the input–output simulations is of the order of seconds on the average laptop.
Performance parameter, 
$\log _{10}({\mathcal{E}})$, as in (6.1), as a function of the weight matrices. (a) IFFC control technique. (b) LQG control technique.
Here, the covariance matrices are constants, as in (3.3), and the weight matrices for the control problem are chosen to be constants,
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \hat{Q}=Q_{IFFC},\quad \hat{R}=R_{IFFC},\\ \unicode[STIX]{x1D64C}=Q_{LQG}\unicode[STIX]{x1D644},\quad \unicode[STIX]{x1D64D}=R_{LQG}\unicode[STIX]{x1D644}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$ In figure 17 the results of the input–output simulations based on the IFFC or LQG methods are shown. From the results based on the IFFC technique it clearly appears that there exists a combination of weights 
$(Q,R)$ where 
${\mathcal{E}}$ is constant. This occurs because the 
$J$ functional in (3.8) can be written as 
$R$, which is a constant in this case, times another functional with only one weight in the form 
$Q/R$. The constant 
$R$ becomes irrelevant in the minimization problem, thus the minimization can be performed with respect to the functional with the weight 
$Q/R$. The weights can be related as
$$\begin{eqnarray}\displaystyle Q=c_{1}R, & & \displaystyle\end{eqnarray}$$ with 
$c_{1}$ a constant. By using (C 3), the solution to the control problem based on the IFFC (3.9) can be written as
$$\begin{eqnarray}\displaystyle & \displaystyle \hat{K}={\displaystyle \frac{{\hat{G}}_{uz}^{H}{\hat{G}}_{yz}}{c_{1}^{-1}+{\hat{G}}_{uz}^{H}{\hat{G}}_{uz}}}. & \displaystyle\end{eqnarray}$$ The results show (figure 17) 
$c_{1}=10^{-4}$. Moreover, by using (C 3), the objective function of the IFFC control technique (3.8) can be written as
$$\begin{eqnarray}\displaystyle & \displaystyle J=R\int _{-\infty }^{\infty }\mathop{\sum }_{k}(c_{1}\hat{u} ^{H}\hat{u} +\hat{z}^{H}\hat{z})\unicode[STIX]{x0394}\unicode[STIX]{x1D6FD}_{k}\,\text{d}\unicode[STIX]{x1D714}, & \displaystyle\end{eqnarray}$$ which shows how under the assumption of constant weights the optimal solution depends only on the ratio of the weights. Writing the cost function as in (C 5) combines the physical meaning of the weights in one single parameter and constrains the solution of the minimization procedure to the isolines 
$J(Q,R)=\text{const}$.
Estimation of TF. Estimating 
$\tilde{\boldsymbol{z}}(t)$ from the available 
$\boldsymbol{y}(t)$. Solid black line, state observer. Dashed red line, empirical TF. (a) Central TF. (b) Complete TF.
Since the LQG results in two independent optimization problems, the control problem and the estimation problem (appendix A), both the objective functions can be expressed in a similar fashion as in (C 5). It follows that the performance 
${\mathcal{E}}$ of the LQG can be expressed as function of two weight matrices only: one weight matrix from the control problem and one weight matrix from the estimation problem. This result is shown in figure 17 as a function of 
$R_{LQG}/Q_{LQG}$ and 
$v_{n}/v_{d}$ (as in (3.3)). Since the variables are associated with two independent optimization problems, there is no general reason for the existence of a set of weights for which the performance parameter 
${\mathcal{E}}$ is constant and has a minimum. In fact, figure 17 shows that 
${\mathcal{E}}$ has a minimum for a specific combination of 
$(R_{LQG}/Q_{LQG},v_{n}/v_{d})$.
The minimum value achieved by the LQG is below the one achieved by the IFFC, 
${\mathcal{E}}_{LQG}^{min}<{\mathcal{E}}_{IFFC}^{min}$. This occurs because the IFFC technique does not include the estimation problem. In figure 18 the estimation function of the LQG corresponding to 
${\mathcal{E}}_{LQG}^{min}$ is compared to the 
$G_{m}^{yz}$, which is used in the IFFC and computed as in § 5.2. It can be seen that the two curves are slightly different. This occurs because the design parameters of the estimation problem can be tuned to seek for the overall optimal solution, which does not appear to be given by the IFFC, even though it was shown that 
$G_{m}^{yz}$ is a good estimation function. In other terms, since the LQG has more design parameters than the IFFC, it can span a space of solutions of higher dimensions than the IFFC. The latter statement also implies that there exists a combination of parameters for the estimation problem in the LQG for which 
${\mathcal{E}}_{LQG}={\mathcal{E}}_{IFFC}^{min}$ holds, and that the result achieved by the IFFC can be seen as a suboptimal solution of the LQG.
The input–output simulations allowed us to identify the weights that give the best performance if the system were to be completely linear.
Open access
