On the wave-cancelling nature of boundary layer flow control

This work deals with the feedforward active control of Tollmien–Schlichting instability waves over incompressible 2D and 3D boundary layers. Through an extensive numerical study, two strategies are evaluated; the optimal linear–quadratic–Gaussian (LQG) controller, designed using the Eigensystem realization algorithm, is compared to a wave-cancellation scheme, which is obtained using the direct inversion of frequency-domain transfer functions of the system. For the evaluated cases, it is shown that LQG leads to a similar control law and presents a comparable performance to the simpler, wave-cancellation scheme, indicating that the former acts via a destructive interference of the incoming wavepacket downstream of actuation. The results allow further insight into the physics behind flow control of convectively unstable flows permitting, for instance, the optimization of the transverse position for actuation. Using concepts of linear stability theory and the derived transfer function, a more efficient actuation for flow control is chosen, leading to similar attenuation of Tollmien–Schlichting waves with only about 10% of the actuation power in the baseline case.


Introduction
The classical route to turbulence for a laminar boundary layer occurs in a low-perturbation scenario, when the mean flow is unstable to infinitesimal perturbations which will grow exponentially in the form of Tollmien-Schlichting (TS) waves. When a critical amplitude is reached, typically of 1% of the unperturbed flow, nonlinear interactions start to occur, eventually leading to three-dimensionality and the breakdown to turbulence, as thoroughly described in the review of Kachanov [27]. The much higher friction drag values of turbulent boundary layers, when compared to the laminar cases, have motivated over the last years several investigations of active closed-loop control targeting at these TS waves.
Nevertheless, flow systems present two characteristics which normally prevent the use of linear control theory to directly tackle the problem; the nonlinearity of the Navier-Stokes equations and the usual high dimensionality, where each grid point of the discretized system may be regarded as a degree of freedom. Such characteristics make the computational requirements both of the simulated system and the online actuation computation to grow rapidly with the dimensions of the calculation domain.
One strategy to deal with these matters consists in designing the control laws off-line, in a reduced-order system, with the ultimate test of the controller being made in the full nonlinear simulation or experiment. This is referred to in the literature as "reduce-then-design" process [44], and has seen several applications, both in numerical and experimental investigations [2,8,42,43]. An interesting reduced-order model for flow control is the Eigensystem realization algorithm (ERA), a technique proposed by Juang and Papa [24], that reproduces the input-output behaviour observed in a simulation or experiment without the need to solve adjoint equations.
Other possibilities include an empirical strategy, using system identification to model the said relationship between inputs and output. Examples of such may be found in the works of Gautier et al. [20] and Hervé et al [22], which make use of ARMAX (Auto-Regressive-Moving-Average-eXogenous), which leads to a time-domain single-input-single-output model, built from unsteady measurements of the fluctuations.
Once the reduced-order model is available, a common strategy for control of boundary layer transition is to place the actuation in a region where the amplitude of the TS waves is small and to account on the convective nature of the flow via a feedforward scheme, where the actuator is placed downstream of the input and upstream of the objective position. The control action is then decided by means of measuring the input and acting to minimize a given quantity at the objective position. This can be accomplished using static compensators, such as the linear-quadratic-Gaussian (LQG) regulator [16,25]. Differently from adaptive schemes, which have been more recently applied due to their increased robustness to plant uncertainties [17,45], LQG is computed off-line, using the reduced-order model, and remains unchanged during its operation on the nonlinear system. Examples of applications of LQG designed using reduced-order models may be found in the works of Barbagallo et al. [6,7] for an oscillator flow and Barbagallo et al. [5] for an amplifier flow. The reader is referred to the work of Schmid and Sipp for an overview of optimal control applied to flow problems [41].
The success of LQG when applied to flow control applications lies on the fact that in spite of relying on an accurate description of the plant, both disturbance and measurement noise can be present and are normally modelled as white noise [4,17]. The states of the system, which are not directly available, are estimated via a Kalman filter; the estimated state is then fed back, and optimal control is obtained by the standard linearquadratic regulator.
In spite of the aforementioned implementations, flow control is yet to be established as a competitive technique and one of the reasons for this fact is due to the energy budget between actuation and reduction in friction drag. With the currently available actuators, the efficiency of closed-loop control is close to the break-even point [16]. Furthermore, the design of the controller based on the ERA model, however efficient, does not provide an interpretation of the states of the system, and therefore, the assessment of the position of sensors and actuators cannot be done without performing simulations using the full system, which, in general, could be computationally demanding.
With this in mind, in this work, we propose a simpler technique, based on the linearity of the Tollmien-Schlichting waves on the region of actuation, to act via a destructive interference of the incoming wavepacket by means of a cancellation signal. Initial works dealing with this sort of wave-cancellation actuation could be tracked down to the first attempts on active flow control; experimental implementations of Milling [35] and Liepmann et al. [33] using a water tunnel with an actuator tuned to cancel a single-frequency wave generated upstream; or Thomas [47], Thomas and Saric [46] and Laurien and Kleiser [29] who, using wind-tunnel experiments and numerical implementations, respectively, were able to act in closed-loop attenuating the TS waves downstream of actuation.
The cited papers, dealt with single-frequency disturbances, with the amplitude of actuation taken from calibration experiments or simulations. The work of Li and Gaster [32] improves on these ideas by formulating a framework in terms of theoretically derived transfer functions (TFs), obtained from the Orr-Sommerfeld equation. The feedforward actuation is again a wave-cancellation signal, and it is obtained by means of an inversion of the estimation and actuation frequency-domain reduced-order models. This accounts directly for the frequency-dependent phase speed and growth rates, and therefore permits the control of broadband incoming fluctuations. Another example may be found in the work of Fan et al. [19] which makes use neural networks to compute the wave-cancellation signal or to model the inverse of the plant, with the implementation of the control law made experimentally.
In this work, we propose to continue to build up on these methods. The theoretically derived, frequencydomain transfer functions are enhanced by using the parabolized stability equations (PSE) solution for this flow, which has the added advantage of accounting for the divergence of the mean flow along the streamwise direction. PSE has been shown over the last years to be a good approach for the prediction of the frequency behaviour of fluctuations in boundary layers, presenting a compelling match in comparison to nonlinear simulations or experiments (see, for example, Bertolotti and Herbert [9], Bertolotti et al. [10] for early works on the use of PSE for boundary layer problems, or Herbert [21], Juniper et al. [26] for an overview of linear PSE for modal analysis). However, PSE had not been used for the development of more control-oriented, time-domain models, until recently [38].
The wave-cancellation control law, derived using the before mentioned transfer functions, is modified by means of an optimization avoiding sensitivities and noise amplifications. The theoretical nature of the wavecancellation scheme derived using the PSE transfer functions, allows the exploration of the parameters of the problem without the need to perform computationally demanding simulations. These properties, in conjunction with the knowledge that the control laws evaluated act by means of a cancellation signal, allow the derivation of a method to determine the most prominent transverse positions for actuation, in an attempt to reduce the actuation cost without degradation of control performance. The ultimate objective of doing so is to go beyond the break-even point, with a significant attenuation of the TS waves and a diminished required actuation.
Finally, the description of the closed-loop problem in terms of estimation and actuation transfer functions permits to evaluate a priori the performance of the controllers in and efficient manner, without the need to implement them in a nonlinear simulations or experiment. This property is explored in this paper to evaluate the penalizations of the kernels and obtain the desired performance.
This paper is organized as follows. In Sect. 2, the configurations used for the simulations are briefly presented. Section 3.1 reviews the concepts of LQG using ERA, and Sect. 3.2 introduces the inversion of the frequency-domain transfer functions. In Sect. 4, a baseline case for the comparison of the control laws is evaluated, along with a description of the method for the choice of penalizations of the controllers and the evaluation of their robustness. Finally, the optimization of the control positions is made in Sect. 5. Conclusions are drawn in Sect. 6.

Flow configuration
The incompressible Navier-Stokes equations model the flow, where u(x, t) = (u 1 (x, t), u 2 (x, t), u 3 (x, t)) and p(x, t) are the velocity and pressure, respectively, at each time step t and position x = (X 1 , X 2 , X 3 ), taken in the Cartesian coordinates. The streamwise direction X 1 is aligned with the free-stream velocity U ∞ , X 2 is normal to the surface of the plate and X 3 defines a spanwise direction. Two-dimensional and three-dimensional simulations are performed; for the former, the dependence of the quantities along the spanwise direction is not considered. A plate of semi-infinite length lies in the X 1 X 3 plane, where no-slip conditions are enforced at X 2 = 0. The control action is analyzed via direct numerical simulations (DNS) with the pseudo-spectral code SIMSON [12], which gives a high numerical accuracy. The same code has been used in several works by this group, both for simulations of laminar/transitional boundary layers [16] and turbulent flows [15].
The flow is periodic along the spanwise direction, which ensures that the boundary layer is essentially two dimensional in terms of statistics, and a fringe forcing, given as λ(X 1 ), introduced in the last 20% of the domain, ensures periodicity also along the streamwise direction. The use of an artificial periodic boundary condition along the axial direction facilitates the use of Fourier methods which, in conjunction with the fringe technique, has been used in the direct simulations of both transitional and turbulent boundary layers (see Nordström et al. [36] for a thorough explanation of the method). The basic idea is to divide the computational domain into one useful and one fringe region, with an extra forcing added to the latter to compensate for the periodicity. Both periodic directions are discretized using a Fourier basis, and for the wall-normal direction Chebyshev polynomials are used. Figure 2 presents a scheme of the simulation and coordinates considered.
Spatial coordinates and velocities are non-dimensionalized using the displacement thickness δ * in the entrance of the domain and the free-stream velocity U ∞ , respectively. The resulting Reynolds number, defined as Re = δ * U ∞ /ν, where ν is the kinematic viscosity, is 1000. Time integration is performed using the fourthorder Crank-Nicolson/Runge-Kutta method, with a constant non-dimensional time step of t = 0.4. The computational domain for the 3D simulation is of [0, 1000] × [0, 30] × [− 70, 70] in the X 1 , X 2 and X 3 directions, with N X 1 = 768 and N X 3 = 96 Fourier modes at the X 1 X 3 plane and N X 2 = 101 Chebyshev polynomials at the vertical direction. The 2D simulation is performed with the same X 1 X 2 configurations.
A volume forcing f is used both to introduce random disturbances and to perform the control action. The forcing is divergent-free, consisting of localized synthetic vortices [43], with spatial support given as indicate the position where the disturbance is centred and χ , γ and ζ the spatial decay of the force. A random broadband signal d l (t) modulates the forcing at the position of each of the upstream perturbations, defined by the index, l (i.e. the perturbations are given by d l (t)b), whereas a control signal u l (t) determines the actuation signal at the position X u , for each of the actuators. Each disturbance is determined independently by uniform white noise, w l (t), and an amplitude, a d , such that d l (t) = a d w l (t). Given that the disturbance forcing is aligned with the stream, with zero component on the spanwise direction, the perturbation will predominantly result in a broadband spectrum of Tollmien-Schlichting waves, with characteristic dispersion relation and group velocity.
In the linear framework, it should be noted that the performance of the controller is independent of the exact amplitudes and frequency content of the perturbations, as long as they result in Tollmien-Schlichting waves at the position of objective. One efficient manner to generate TS waves is by using the chosen spatial support, given in Eq. 3. Furthermore, the impulse response of the perturbations will later be used in an intermediary step, in the construction of the transfer functions and Eigensystem realization algorithm matrices, which require the impulse response of the perturbation present on the system. A localized measurement of the streamwise skin friction is used to define the inputs given by sensors y(t), and downstream objective, z(t). Four rows of equispaced independent objects are placed with a transverse separation of X 3 = 10 for the 3D simulation; Semeraro et al. [43] show that such spanwise spacing is necessary to the current setup, given the transverse wavenumber content of the output. Disturbances are added at the wall, at X 1 d = 35, and measurements are taken at X 1 y = 300 and X 1 z = 550, defining input and objective, respectively. Actuation is performed at X 1 u = 400, both on the wall, for the baseline case, an at an optimized position, derived in Sect. 5. A simplified scheme of the set is shown in Fig. 3.
The axial separation between the actuation and objective was determined by evaluating the spatial transient induced by the actuator, as will be shown in further details in Sect. 5.

Control methods
Through this section, the two control strategies that were evaluated will be introduced. Both of them are based on a feedforward scheme, which is in accordance with the convective nature of this flow, where downstream actuation u has negligible effect in the upstream sensor y [8]. The resulting block diagram, written with the variables transformed to the frequency domain, is shown in Fig. 3 and makes it clear that the resulting output is given by a superposition of an open-loop behaviour, determined by the transfer function G yz , with the actuation, determined by the transfer function G uz . The strategy to calculate the convolution kernel K will define the two evaluated control laws: linear-quadratic Gaussian (LQG), which is designed using the state-space representation of the system, and inversion-based control, which deals directly with the frequencydomain transfer functions of the problem.

LQG control
LQG controllers have been successfully applied in flow control problems over the last years. For applications in boundary layer transition control, the reader is referred to the works of Fabbiane et al. [16,18], and the work of Bagheri et al. [2] which presents a nice introduction to optimal control for applications in fluid mechanics. The success of LQG controllers is related to their optimality, the design being made in two steps, from the solution of two Riccati equations for the estimation and control problems, which guarantee stability, as long as the system is observable and controllable.
The control scheme is based on the theory of optimal linear-quadratic regulators (LQR). The objective of the controller is to find an actuation signal u(t) in order to minimize the energy of a given state, q(t). The energy budget of the actuation is also considered, such that a cost functional of the form of Eq. 4 results.
Q and R are penalizing matrices for the state and actuation, respectively, and the subscript H indicates the conjugate transpose of a given matrix. The minimization of the cost function is subject to a restriction given in terms of the state equation, where A and B represent the system dynamics and the effect of an actuation. The Lagrangian is then defined as considering the functional and the restriction where the adjoint state h H plays the role of the Lagrangian multiplier. The condition for optimality is given by the derivative of the Lagrangian with respect to the actuation, as per equation 7, such that is the optimal control signal, which minimizes the previously defined cost function, subject to the restriction given in terms of the state equation. A linear relation is assumed between the direct and adjoint state, h(t) = X(t)q(t), such that the gain results in The matrix X(t) is a solution of a Riccati differential equation (see Lewis et al. [30]), and it is, therefore, possible to calculate it off-line and only once, prior to its implementation on the nonlinear problem. Solvers for the Riccati equation are available in standard software for low-dimensional problems [1]. This implies that a model reduction should normally be performed for applications in flow control problems. Furthermore, in flow control applications, the full state q(t) is usually not available, and therefore has to be estimated from an finite number of available measurements via an estimation problem which in LQG is addressed via a Kalman filter. We start by considering the state-space representation of the observer, as per Eq. 10.q = Aq + Bu(t) − L(y(t) −ŷ(t)) whereq(t),ŷ(t) andẑ(t) are the estimations of the state, measured input and output, given in terms of the matrices C y and C z . The objective is then to design the filter L minimizing the error between the real and estimated state subjected to stochastic white noise with zero mean in the form of disturbances and measurement noise. The cost function for minimization is based directly in the error. The filter L is also obtained through the solution of a Riccati equation [2]. The control law can then be written combining the Kalman filter and LQR, feeding back the estimated state [4], with the input as the measured signal y(t) and the output the control response, u(t).
leading to a convolution kernel in the form of a finite impulse response filter (FIR), which, in the discretized version reads where N determines the buffer for calculation of the convolution. Extension to the 3D problem is made via a multiple-input-multiple-output (MIMO) strategy [16]. Given that the control law is linear and there is an equal number of sensors and actuators, (M + 1), the total number of transfer functions to be considered, in the 3D case, would be (M + 1) 2 . However, the flow is homogeneous and periodic along the spanwise direction, which implies that the same set of transfer functions between one actuator and the respective sensors can be replicated for all the others [16]. Such assumption reduces the number of transfer functions to (M + 1), and the control law for each actuator is calculated by a double, discrete convolution, as defined in Eq. 14 where M + 1 is the number of sensors. LQG is found to be very effective when a precise description of the system is available, whereas the exact form of the perturbations that are present is not exactly known. The fact that the resulting actuation is optimal in terms of the budget between actuation and objective is also a strong suit of this method. Nevertheless, the exact mechanism which results from the actuation is not immediately clear, particularly as for this case, the flow physics are replaced by the ERA matrices, which cause a physical interpretation of the states to be no longer available. This will be further discussed in the next section

Reduced-order model: the Eigensystem realization algorithm
To obtain the LQG controller, in general two Riccati equations should be solved; the first to obtain the Kalman filter, and the second to obtain the control law [4]. This is only feasible if the order of the system is not too high. For fluid flows, a reduction of order is, thus, required [3,28]. In particular, a reduced-order model should be able to reproduce the input-output behaviour of the flow, with a good representation of the transfer functions from inputs d and u to outputs y and z. A schematic representation of such transfer functions is shown in Fig. 3.
To compute the controller, the influence of the perturbation and actuation at X 1 d and X 1 u into the inputs and outputs of the system has to be available in the form of a state-space relation. This can be done via the Eigensystem realization algorithm, which was initially proposed in the work of Juan and Pappa [24], and leads to a model that reproduces the impulse responses of the system with reasonable accuracy. The algorithm is described in detail by Ma et al. [34], and also applied to the flow over a flat plate at high angle of attack. The method is based on a set of impulse responses assembled to build a Hankel matrix, whose singular value decomposition supplies the state-space matrices which reproduce reasonably accurately the impulse responses of the real problem. For completeness, a brief description of this method is outlined in "Appendix".
We obtain ERA models for the boundary layers considered here by simulating impulse responses in time, with the defined spatial support of the forcing, from the perturbation and actuator. For illustration purposes, the impulse responses of the 2D simulation in comparison with the resultant ERA model are shown in Fig. 4. The resulting model has N ≤ 100 degrees of freedom, which is computationally treatable.
As shown in Fig. 4, ERA system identification reproduces accurately the dynamics of the system. A strong characteristic of this method is that it does not require the use of the adjoint system and is, therefore, appropriate both to numerical simulations and experiments. It has recently been shown [34] that there is a theoretical equivalence between ERA and balanced proper orthogonal decomposition (POD). This implies that the former produces models that are approximately balanced, truncating states whose controllability and observability are lower than a given established tolerance. Even though this property supplies a physical interpretation to the resulting reduced-order model, at least for what concerns problems with a high number of degrees of freedom, such as boundary layers, there is no direct interpretation of the states of the ERA matrices in terms of the original (not reduced) physical problem.

Control law
For this convectively unstable flow, with the actuator downstream of the input, in a feedforward configuration, as shown in Fig. 3, a straightforward strategy for control is to perform the direct superposition of the open-loop wavepacket with the actuation, aiming at a direct cancellation of the resulting output. This is made in the frequency/spanwise wavenumber (ω/β) domain 1 , as where Z ol (ω, β) and G uz (ω, β)U (ω, β) are the open-loop behaviour and the effect of actuation, respectively, written in terms of their ω/β form. The transfer functions G yz (ω, β) and G uz (ω, β) are, respectively, the estimation and actuation transfer functions and account for the open-loop behaviour of the system and the effect of an actuation in the output position. G yz relates two outputs of the system, the sensors y and z. Such transfer function can be obtained if the system is dominated by a single amplified mode travelling between the sensors [38]; this is the case for the TS-dominated transition studied here.
The frequency-dependent phase speeds and growth rates of the Tollmien-Schlichting waves are readily considered in these transfer functions, and therefore, the kernel also has such dependence. We, thus, have U (ω, β) = K (ω, β)Y (ω, β) in the transformed plan, and the resulting output may be written as The gain may be calculated via the direct inversion of the actuation and estimation transfer functions to obtain a null output, such that, The resulting kernel, in the time domain, is obtained from the inverse Fourier transform (from ω to t and β to z) of Eq. 17, leading to the same actuation signal as in Eq. 14. This procedure has been successfully applied to boundary layer control in the work of Li and Gaster [32], where the transfer functions were derived from linear stability theory. The inversion has the drawback that it will be ill-conditioned in the zeros of G uz (ω, β), which may be interpreted as uncontrollable frequencies and wavenumbers of the problem, for which a Tollmien-Schlichting wave is stable at the position of actuation. Furthermore, this method does not account for plant uncertainties or penalize for high actuation values, which may result at such uncontrollable cases. This problems have been tackled in the work of Devasia [14], with a systems theory perspective. The basic idea is to deal with a functional defined in the frequency domain, where R(ω, β) and Q(ω, β) are penalizing matrices. Due to Parseval's theorem, this method is analogous to the linear-quadratic regulator (LQR), but defined in the frequency, rather than in the time domain; a difference is the capability of setting ω and β dependent penalization terms. It should be noted, however, that a direct comparison with the LQG method is not possible as an observer is included in the system by means of a Kalman filter for this case. Minimization of Eq. 19 leads to the penalized inverted gain [14], When compared to the LQG framework, the advantage here is that there is a clear interpretation of the control action. It results in a cancellation signal, leading to a destructive interference between the incoming TS wavepacket and the TS waves excited by the actuator. The states of the problem are directly obtained from the frequency-domain transfer function, such that there is no need to add an extra step given by the Kalman filter estimation. The requirement to apply this technique is the availability of the transfer functions in Eq. 20, and the methods for which will be described next.

Calculation of transfer functions
Given that the wave-cancellation gain is posed uniquely in terms of frequency-domain transfer functions, it is possible to deal directly with the impulse responses of the perturbation and actuator, written in their (ω/β) counterparts. This avoids the use of ERA matrices. The TF from u to z is obtained directly from the Fourier transform of the impulse signal taken from a linear simulation, whereas the TF from y to z is obtained by the quotient of the TFs from d to z and d to y, respectively, as given in the equation below, The perturbation impulse response here is used on an intermediary step in the determination of the estimation transfer function, G yz (ω, β), which is insensitive to the exact details of the perturbation, as long as it is located upstream of the two measurements and leads to Tollmien-Schlichting waves at their position. To illustrate this method, the impulse responses from the actuators, along with the constructed TF, from y to z are shown in Fig. 5, for the 3D problem. Analogous results may be readily obtained for the simpler, 2D case. The impulse responses in Fig. 5 lead to a Gaster wavepacket, related to the fact that each β/ω mode has a different growth rate, associated with the corresponding dispersion relation. The time delay corresponds to the distance between the two positions divided by the group velocity, which is of approximately 0.4, for this flow.
The frequency content of the impulses presents a maximum at β = 0 and at the TS most amplified frequency for this region of the flow. The region where the G uz(ω,β) tends to zero corresponds to the uncontrollable frequencies of the problem and has to be either filtered or removed via the optimization of the functional in Eq. 19. It should be noted that the β/ω content of the actuation should lie in approximately the same region of that corresponding to G yz (ω, β) such that the system is controllable and avoiding high actuation values.
Since the only necessary ingredients for implementation of wave-cancellation control are the estimation and control TFs, an alternative to obtain these is to proceed theoretically and obtain the transfer functions using the linear parabolized stability equations (PSE) [21], initialized with the solution of a locally parallel problem. This idea was firstly introduced in the work of Sasaki et al. [38] with the objective to estimate the pressure fluctuations in a high-Reynolds turbulent jet, and it was used as a reduced-order model for active control of a mixing layer in [37,39]. For the sake of clarity, the method will be briefly outlined here; for a more complete description of the method as a time-domain predictor, the reader is referred to one of the above cited works.
The PSE solver considers a pseudo-spectral method used for the discretization of the functions along the wall-normal direction by means of Chebyshev polynomials. The axial derivatives are calculated using finite differences, and the spatial marching is performed with an implicit Euler method. A mapping is applied to transform the Chebyshev domain to an infinite one, where Dirichlet boundary conditions are imposed at X 2 → ∞. The system resulting from this operations may be found in Juniper et al. [26]. For details concerning the nature of the PSE solution or a review of its use for modal analysis of boundary layer flows, the reader is referred to the works of Li and Malik [31], Herbert [21] and Juniper et al. [26].
In order to construct the models using PSE, the idea lies in finding the estimation and actuation transfer functions, G yz (ω, β) and G uz (ω, β) by solving the PSE for wide range of frequencies. This is feasible due to Time-and frequency-domain transfer functions in the first and second columns, a, b, d are actuation TFs extracted directly from the linear simulation, from d to y and z and from u to z. c is the constructed estimation TF from y to z. a Impulse response from d to y. b Frequency/transverse wavenumber content of the impulse from d to y. c Impulse response from d to z. d Frequency/transverse wavenumber content of the impulse from d to z. e Constructed impulse response from y to z. f Frequency/transverse wavenumber content of the impulse from y to z. g Impulse response from u to z. h Frequency/transverse wavenumber content of the transfer function from u to z the low computational cost involved in the method. For the estimation case, we have, where the quantities Z (ω, β) and Y (ω, β) are taken from the PSE solution, and the quotient eliminates the free constant inherent to the solution of this linear method. To obtain the TF in the z and t variables, the inverse Fourier transform (analogous to Eq. 18) of G yz has to be taken. The procedure to obtain the actuation transfer function, G uz , follows a somewhat different procedure, which is outlined in detail in [40] and consists in starting the PSE upstream of the actuation, such that transients are able to be damped. The value of the envelope of the solution at the position of actuation is then corrected using the projection of the forcing into the Tollmien-Schlichting mode. This is performed for a wide range of (ω, β), Fig. 6 Comparison between the prediction obtained from the PSE TF and a nonlinear simulation of a 2D boundary layer and the method to obtain the time-domain counterpart of the transfer function follows exactly the same ideas outlined for the estimation problem.
After having the TFs available, the computation of the control law follows the same path as when the impulse responses are used as reduced-order models. Given that PSE is capable of following a single mode, it is started with a Tollmien-Schlichting wave and only perturbations associated with it will be predicted by this method. As an example, the PSE prediction is shown in comparison to a nonlinear simulation of a 2D boundary layer in Fig. 6, for the output at X 1 z = 500 and the input at X 1 y = 300, with broadband perturbations inserted at X 1 d = 35. The PSE prediction is obtained from a convolution of the measurements y(t) with g yz (t), A compelling agreement between the reduced-order model and the simulation is observed. The use of a theoretical method such as PSE to compute the control law permits the exploration of certain parameters, such as the optimal transverse position for the actuators, prior to the computation of the kernel and without the need to perform computationally demanding simulations for each configuration. The method to determine the most efficient transverse position will be described in detail in Sect. 5.

Baseline case
The baseline case for study considers the positions of 300, 400 and 550 for input, actuation and output (X 1 y , X 1 u , and X 1 z ), respectively. The separation of 150 between actuation and output was chosen such that the latter was downstream of the transient induced by the actuator, both when acting at the wall (meaning that X 2 0 = 0, in Eq. 3, and the resulting spatial support of the forcing is such that its peak is located at X 2 = 0) or at the optimized transverse position (X 2 = 0), which will be derived in Sect. 5. Figure 7 shows the resulting RMS value of the axial velocity fluctuations shear at wall when introducing a broadband forcing the position of X 1 u = 400 at the wall (X 2 0 = 0) and at the optimized transverse position. It should be noted that the actuation is introduced in open loop, without any other perturbations (d(t) = 0). Its initial fast growth is related to a spatial transient. Once the transients cease, the velocity fluctuations start to behave exponentially, on the characteristic trend of a Tollmien-Schlichting wavepacket, which then results on a linear growth when evaluated in a logarithmic plot.
It is noticeable that the spatial transient, during which the fluctuations are not behaving exponentially, induced by the optimized position is longer, the position of 550 being already in the region of exponential growth for both cases, where Tollmien-Schlichting waves, which will be targeted by the control law, dominate the dynamics.

2D simulation
We start by obtaining control laws for the two-dimensional simulation. Prior to the implementation of control, it is necessary to determine the parameters for the minimization of the functionals. The choice of penalization for LQG and wave-cancellation was made by fixing the penalization of the objective at unity and varying the penalization of the actuation in order to obtain a regularization of the resulting kernel and the highest attenuation at the output. The process is analogous between LQG and wave cancellation, both in 2D and 3D, and it will be thoroughly explained in Sect. 4.2.1, for the more complicated, 3D case. The resulting values of this method were R = 1 and R = 100 for LQG and wave cancellation, respectively. The wave-cancellation penalization was taken as constant for all frequencies and wavenumbers. As for the LQG penalization, it was considered the same for the different states.
Wave-cancellation kernels were calculated using the PSE and impulse responses transfer functions. Figure 8 shows the behaviour of the resulting kernels in the time domain along with their frequency content, given in terms of the power spectral density. Differences may be seen between the three kernels mostly in the frequency domain, and in the low frequencies, where LQG is disposing of less energy than the other two strategies. This trend is most likely to occur as the actuator is unable to efficiently act on these frequencies, and they are less present on the output than the higher ones, more characteristic of TS waves. Figure 9 shows the behaviour of three closed-loop cases, when compared to the uncontrolled simulation, in time and frequency domains. A compelling reduction of more than one order of magnitude is observed for all the kernels, with LQG significantly out-performing the other two. This should be understood due to its kernel being calculated from an optimality condition taking into account both the energy of the output and actuation. The ERA reduced-order model also leads to an accuracy which is as high as the one obtained from the TFs using the impulse responses of the system. Differences concerning the two wave-cancellation techniques only occur due to the PSE estimation and actuation transfer functions leading to small inaccuracies in the prediction of the output, when compared to the direct use of the impulse responses of the system. These inaccuracies lead to small differences in the kernels (Fig. 8), particularly on what concerns their phase. Such differences account for the different performances observed in the closed-loop system (Fig. 9), with the inversion based on the impulse responses leading to more compelling reductions.
There is a strong similarity between the wave-cancellation and LQG kernels, particularly when the impulse responses are used for the computation of the former, which shines a light on the mechanisms of control of LQG indicating that they approach those of wave cancellation. A compelling amplitude reduction of more than two orders of magnitude in the objective sensor z is obtained for the most amplified frequencies. A similar reduction is also observed in the root-mean-square (RMS) value of the fluctuations.

Choice of penalizations
This section has the main purpose of demonstrating how the choice of penalizations for the optimal controllers is made. Since the process is the same in 2D and 3D cases and performed analogously for both LQG and wavecancellation (WC) kernels, we have chosen to demonstrate the method for the wave-cancellation kernel applied in the 3D problem, only. Furthermore, WC kernels have an extra characteristic which makes the choice of penalizations more critic: the inversion (Eq. 17) could lead to a noisy kernel in ω, β regions where the actuation transfer function, G uz (ω, β), presents a low amplitude. These correspond to uncontrollable frequencies of the problem, and the resulting kernel should not have a significant gain at them in order to avoid the amplification of unwanted noise. The inclusion of the penalization term in the determinations of the WC kernels, shown in Eq. 19 alleviates this problem, provided Q and R are chosen in an appropriate way. The evaluated controllers present two types of penalizations: in the objective state and on the actuation signal, Q and R, respectively. Penalizations are then sought such that the resulting kernel is regular, without spurious behaviour, and also lead to the desired performance in closed-loop. The resulting sensitivity of the closed-loop system is then verified a posteriori.
It was chosen to evaluate these characteristics via two methods. In order to evaluate the regularization, we define following Cordier et al. [13], where ±β n is the Nyquist frequency, based on the discretization along the transverse direction. The parameter L is a measure of the amplitude of the kernel and should be evaluated for several values of either Q or R. Low values of L are expected if control penalization is high, leading to low actuation gains. On the other hand, insufficient control penalization might lead to high values of L, with the actuation occurring for frequencies and wavenumbers with low controllability. It was chosen to keep the value of Q at unity and vary the value of R. The second metric considers the expected performance of the kernel, when applied in the closed-loop problem. This test may be done directly in the system described in terms of transfer functions, the inverse Fourier transform of Eq. 16, Estimating the output z in the closed-loop system based on Eq. 25 is computationally inexpensive when compared to the test in the direct numerical simulation, and allows the evaluation of several kernels and therefore a choice of penalization that would have the desired performance in the closed-loop nonlinear problem. Figure 10 shows the behaviour of parameter L and the closed-loop performance of the resulting kernel (tested using the transfer function-based framework), for different values of R. The different regions highlighted in Fig. 10 correspond to considerably different kernels, and examples of such are presented in Fig. 11.
In I , which corresponds to low values of penalization R, the resulting kernel has high and fast changing amplitudes, due to the fact that the control is acting at frequencies for which the actuator is unable to efficiently generate flow fluctuations; moreover, such frequencies are not significantly present in the output. Application of this kernel in closed loop leads to noise amplification and therefore to an increase in the output, rather than a decrease. By continuing to increase the penalization, the L parameter tends to stabilize, between regions I I and I I I , with much smoother variations, where it starts to lead to a decrease in the output. However, a high effectiveness has not yet been achieved. In region I V , the most balanced kernels are reached, which are able to effectively cancel the incoming TS waves, without amplifying undesirable frequencies. Only in region I V , we start to observe kernels with peak energies at the (ω, β) combinations corresponding to amplified TS waves, as can be seen by the crosses in Fig. 11. If one continues to increase the penalization in the objective, the effect is mostly on the amplitude of the resulting kernel. This corresponds to region V , with a corresponding closed loop without amplification of noise, but which also does not attenuate the incoming TS waves efficiently.
Therefore, a kernel in region I V was then chosen to be implemented in the nonlinear closed-loop problem. The same process was applied to the LQG design, keeping the penalization on the objective constant and varying that of the objective, while evaluating simultaneously the regularization of the kernel and its performance in closed loop.

Closed-loop results
The analysis is now performed in the nonlinear 3D simulation described in Sect. 2. The resulting wavecancellation and LQG kernels account for the 15 sensor/actuator pairs and are shown in Fig. 12. Both the continuous and discretized versions of the kernel are shown in this figure. The continuous version of the kernels would permit the calculation of the actuation via a 2D convolution with the input measurements u(X 3 , t) = ∞ 0 The kernels are implemented in the simulation using their discretized format, which permits the direct use of Eq. 14.
Analogous to the 2D case, the similarity of the two kernels obtained by the different methods is compelling, supplying further evidence that optimal control obtained via LQG is also acting by means of a wave-cancelling signal.
The behaviour of fluctuating wall shear stress, with respect to the unperturbed, laminar solution, is shown in Fig. 13 for the uncontrolled and controlled cases. The positions of perturbation, input, actuator and objective are the same as for the 2D case. As expected, the perturbations grow exponentially while being convected from their generated positions, at X 1 d , into a shape typical of a Tollmien-Schlichting wavepacket.
There is a clear reduction in the amplitude of the TS waves downstream of the objective position, for both the evaluated strategies. The resulting actuation of the two controllers appears to be equivalent, with the same type of structure and amplitude remaining downstream of the actuator.
The ω/β content of the three cases, given in terms of their Fourier transforms, is shown in Fig. 14. The peak at the uncontrolled case is located at β = 0, at a frequency correspondent of TS waves, confirming that these dominate the flow for the choice of perturbation d(t). This peak is significantly attenuated for both controllers, with almost no differences in terms of their performance. Compelling reductions of about 75% of the objective signal were observed for the use of both controllers. The same trend was achieved in the most amplified frequencies and transverse wavenumbers. A similar performance was again observed between LQG and wave-cancellation kernels, as was shown for the 2D cases.   The robustness of the two derived feedforward controllers was evaluated with respect to variations of three parameters: amplitude of inflow perturbations, a d , the presence of noise in the input sensor, y(t), or in the closed-loop actuation, u(t). The results shown in this section will concern the more challenging, 3D case; similar results are expected to be seen in the simpler, 2D case. Increasing the amplitude of the input perturbation tends to trigger nonlinear interactions sooner, which in turn is expected to degrade the performance of the controllers, since the derivation follows a linear assumption. Figure 15 shows the variation of the quotient between the mean square values of the controlled and uncontrolled signals, at the positions of objective, X 1 z , when the amplitude of perturbation is increased. It is noticeable that when a d becomes higher than 3 × 10 −3 , the performance of the two evaluated control laws decreases rapidly.
Such tendency may be understood in more depth from Fig. 16, which shows the amplitude of the measured signal y(t), at the input position, X y , for an increasing value of a d . It may be observed that for values of a d higher than 2 × 10 −3 , the measured signal no longer presents a linear scaling with a d , indicating that nonlinear behaviour starts to play a role in the dynamics of the system. This trend permits to understand the resulting behaviour of the controllers.
The other two evaluations concern the inclusion of noise in the measurement y(t) and actuation u(t), and seek to model a more realistic experimental implementation, where environmental noise could be present For both cases, the resulting actuation and measurement were taken as a superposition of their real values with a white noise perturbation of fixed amplitude, Eqs. 27 and 28.
where u real and y real represent the noiseless values of the parameters. Figure 17 shows the behaviour of the controlled and uncontrolled cases for a varying amplitude of measurement noise until 50% the amplitude of its unperturbed value. A similar behaviour is seen for both the controllers when noise is added to the actuation and will not be presented, for the sake of brevity. It is noticeable that both controllers are fairly insensitive to the addition of noise, with LQG slightly outperforming the wave-cancellation controller, which is expected given that its design already considers the presence of measurement and actuation noise. It should also be observed that the addition of white noise is only expected to affect the closed-loop performance of the controller in the narrow band corresponding to the unstable frequencies of the Tollmien-Schlichting waves.

Improvement of the transverse position for actuation
The results of previous sections indicate that the closed-loop control is acting by means of a wave-cancelling signal of the incoming TS wavepacket, such that destructive interference occurs downstream of the actuator position. With knowledge of this fact, it is now possible to use the PSE TFs and concepts of linear stability theory to determine prominent positions for actuation in order to increase the efficiency of the controller. This will be evaluated in the simpler 2D case.
The idea is to use the position of the critical layer, where the receptivity of the flow is greatest to perturbations [23] and then consider the different amplification rates of each frequency by means of the PSE transfer function. By doing so, a position is chosen such that TS waves are efficiently generated in the actuation position, and are mostly amplified at the output region.
The transverse position of the critical layer, X 2 crit , varies for each frequency and axial position and is located where the phase speed of the fluctuation is equal to the velocity of the unperturbed base flow. The position of the critical layer, as a function of the angular frequency of the perturbation, at the position of actuation, is shown in Fig. 18. Acting at these positions is expected to be effective in generating TS wavepackets at the corresponding frequencies.
However, each frequency is amplified at a different rate at the position of objective and the precise determination of the optimal position for actuation, X 2 opt , should include this fact. We start by defining the parameter in Eq. 29, which is a function of the transverse position of actuation, Y , and corresponds to the integral of the PSE transfer function over the unstable frequency range, ±ω uns , the parameter (X 2 ) accounts for the different amplification rates of each frequency and will be used to perform an average weighting of the position of the critical layer. The resulting optimal transverse position for actuation is shown in Eq. 30. Figure 19 shows the behaviour of parameter , with the corresponding improved position for actuation, X 2 opt , obtained from Eq. 30, highlighted by the transverse line. The behaviour of the controller is going to be evaluated at this position. It should be noted that perturbations close to the critical layer lead to longer transients, as shown in Fig. 7; therefore, all the comparisons were made at the objective position of X 3 = 550, which is downstream of the transient for the two transverse positions of actuation considered.
The wave-cancellation and LQG kernels were computed for the same objective position but considering the optimized transverse position of actuation. To illustrate the effect of the change in the transverse position in the kernel, Fig. 20 shows the time-domain wave-cancellation kernels for actuation on the wall and at the optimal position. We observe a clear reduction in the amplitude of the kernel for actuation at X 2 opt , which is expected to lead to lower actuation values, since the input position was kept constant.  In order to confirm the reduction in the resulting actuation values when acting at the optimized position, LQG and wave-cancellation kernels were computed and evaluated on the 2D DNS. Different penalization matrices (Q and R) were used for the LQG case and chosen such that the amplitude of the resulting output was approximately the same, for the two actuation positions evaluated. Table 1 summarizes the findings for  the different cases.  Table 1 indicates there is a reduction in the actuation power close to two orders of magnitude, without a loss of performance of the LQG controller, when acting at the improved transverse position, determined from Eq. 30. Figure 21 presents the rms of the axial velocity fluctuation as a function of the axial position and confirms the reduction in the axial velocity remains downstream of actuation. The same behaviour is also observed for the wave-cancellation control.
The changes made to the actuator, in terms of its spatial support, are significant-the peak of the forcing was moved to a position above the wall, chosen in accordance with the method outlined in this section. The design of the kernel for closed-loop control accounting for this different forcing allowed a similar performance with a much lower actuation power. The results of this section tackle the limitations of transition control highlighted by Fabbiane et al. [16], who showed that plasma actuators dispose of roughly the same energy as that saved by the transition delay. The improved actuator is capable of efficient Tollmien-Schlichting wave generation, leading to a cancellation with reduced required actuation power, which improves the overall performance.
Nonetheless, it should be noted that with the currently available actuators for closed-loop flow control it is not possible to act with a spatial support that has a peak above the wall, particularly when dealing with a streamwise actuation, which is necessary to the effective excitation of TS waves (see Cattafesta and Sheplak [11] for an overview about actuators for active control). This section, however, serves as a proof of concept, showing that there would be a way to obtain compelling gains in performance by changing the spatial support of the forcing along the normal direction, trying to obtain the highest possible value at the most sensitive position. The outlined method could also be used as a guideline for the design of plasma actuators or aid in the choice of actuator most suitable for a given implementation. The wave-cancellation approach and the framework described in terms of transfer functions are fast methods to evaluate the potential benefits and different performances of new actuators, without the need to execute computationally demanding simulations or experiments.

Conclusions
In this work, we introduced a technique, building upon the work of Li and Gaster [32], and derived in the frequency domain for the closed-loop feedforward control of fluctuations in 2D and 3D boundary layers. The control law, which leads to a wave-cancellation signal, was modified by means of and optimization, performed in the frequency domain, and firstly described in Devasia [14], which permitted to avoid high sensitivities and noise amplifications. Two reduced-order models were considered as necessary for the method; impulse responses obtained from a linear simulation and PSE transfer functions. The more theoretical PSE TFs show a reasonable agreement with nonlinear simulations and the Fourier transform of the impulse responses.
Results were compared to the standard linear-quadratic-Gaussian regulator (LQG), designed for a reducedorder state-space model obtained through the Eigensystem realization algorithm (ERA), and similar control performances were obtained. Reductions of two orders of magnitude in the RMS values of the shear induced by the axial velocity fluctuations were observed in the 2D problem, and attenuations of 75% were obtained in the nonlinear 3D case.
The evaluation of the robustness of the controllers indicated similar performances; both controllers are able to efficiently control the flow, leading to compelling reductions of the objective, as long the amplitude of the seeded perturbation does not trigger nonlinear effects. When evaluated in a noisy environment, via the addition of white noise to measurement and actuation, both controllers behaved well, leading to an effective noise reduction, indicating the possibility of their implementation in real experimental implementations.
Each strategy presents different interesting characteristics which are relevant for control. The simplicity of the inversion kernels permits a direct interpretation of the actuation signal, and their derivation using the impulse responses avoids the introduction of the ERA system, whose state no longer has a physical interpretation. On the other hand, the use of the theoretically derived PSE transfer functions, although slightly less precise and consequently with a lower performance, allows the derivation of kernels without using computationally demanding simulations. Exploration of the parameters involved in the control, such as the transverse position of actuation, is also possible using this method. The LQG regulator, as expected, given its optimality, presented the best results in terms of attenuation of the objective.
The physical interpretation of the wave-cancellation kernels, and the possibility to derive them from the theoretical PSE transfer functions, enabled the possibility to find an improved position for actuation. This was performed by obtaining the position of the critical layer, where TS waves are most sensitive to forcing, combined with the consideration of the different amplification rates for each frequency, as predicted by the PSE transfer functions. A compelling reduction in the actuation signal was obtained without a significant loss of control performance.
Finally, the similarity of the evaluated controllers allows insight into the physics behind the active control of boundary layers, indicating that the control signal acts by introducing a wave-cancellation signal which destructively interferes with the incoming TS wavepacket, leading to considerable reductions in the fluctuations downstream of the objective.