Extremum seeking to control the amplitude and frequency of a pulsed jet for bluff body drag reduction

Feedback control of fluid flows presents a challenging problem due to nonlinear dynamics and unknown optimal operating conditions. Extremum seeking control presents a suitable method for many flow control situations but involves its own challenges. In this paper, we provide a brief analysis of the extremum seeking method, with attention to modifications that we find to be advantageous. In particular, we present an adaptation for optimisation of the frequency of a harmonic input signal, a common scenario for open-loop flow control systems. We then present results from the experimental implementation of our modified method to the open-loop control system of Oxlade et al. (J Fluid Mech 770:305–318, 2015), an axisymmetric bluff-body wake, forced by a pulsed jet. We find that the system is able to achieve optimal operating conditions in both the amplitude and frequency of the harmonic input signal, and is able to largely reject the disturbances arising from measurements of a highly turbulent flow. We finally show the ability of the extremum seeking system to adapt to changing conditions.


Introduction
Fluid flows are ubiquitous throughout systems of engineering interest, and the behaviour of such flows is often of critical importance. For example, in heat exchangers, the properties of the flow determine the heat transfer, while in the flow over a bluff body, the flow determines 1 3 159 Page 2 of 14 to changing conditions. In such situations, it is often desirable to simply seek to maximise or minimise some parameter such as drag, mixing or heat transfer, adapting to a new optimum as conditions change. In such cases, an extremum or slope-seeking controller may be appropriate.
ES control is a form of adaptive closed-loop control that acts as a real-time optimisation method. This is appropriate for systems with a nonlinear equilibrium map, valid in a time-averaged sense, on which an unknown target location may be identified by its gradient. In this case, no model is required for the flow, only the assumption of a peak or location of known gradient in the equilibrium mapping between a reference parameter and the output. For example, such a mapping could be that between the angle of attack and the lift of an aerofoil, which exhibits a peak at a certain attack angle. ES control was first suggested in 1922 (see Tan et al. 2010, for a review), but has in the last decade or so seen a resurgence of interest following proofs of stability by Krstić and Wang (2000) and Krstić (2000). Included in this recent interest has been the application to a number of fluid flow systems. For example, ES control has been used to directly control the parameters involved in open-loop control systems: Becker et al. (2007) and Pastoor et al. (2008) controlled the amplitude of pulsed jets used on an aerofoil and bluff body, respectively. Alternatively, ES control can be used to adapt some parameter within a closed-loop system. For example, Kim et al. (2009) used an ES controller to adapt the phase shift of the feedback gain used in the control of a cavity flow.
While a standard algorithm exists for ES, many previous authors have suggested modifications, including in the context of flow control. For example, Beaudoin et al. (2006) used a variant on the traditional gradient estimation, computing a real-time fast Fourier transform as an alternative to the often used demodulation (see Sect. 2.1). Alternative gradient estimation methods have also been suggested by Henning et al. (2008) who used an extended Kalman filter, finding it to be much faster than the estimation of the standard algorithm. This algorithm was also tested by Wu et al. (2015) who confirmed this result. Authors have alternatively looked at methods to adapt the control gains. For example, Chabert et al. (2014Chabert et al. ( , 2016) used a fuzzy-logic method to adapt the gains according to the current conditions. In general, many different modifications may be advantageous; however, their suitability depends upon the particular system to be controlled.
While the application of ES to fluid flows has been demonstrated to be effective, there are considerations pertinent to flow control which have not previously been addressed. Principle among these is the challenge involved in optimising the frequency of a harmonic signal. In open-loop flow control, this is a particularly common scenario. For example, in the use of synthetic jets (Kourta and Leclerc 2013;Kotapati et al. 2010), pulsed blowing (Seifert et al. 2008;Chabert et al. 2016), oscillating tabs (Bigger et al. 2009) or in a number of other situations (Greenblatt and Wygnanski 2000), the control consists of a periodic forcing. In many such cases, performance of the open-loop system is dependent on the frequency of the periodic signal, the optimal frequency depending upon the particular control mechanism and the particular operating conditions. There are therefore advantages to be gained through the application of ES controllers, both to find an optimal frequency and to adapt as conditions change. We will show that frequency optimisation of a harmonic signal cannot be done with a standard algorithm and propose a modified algorithm suitable for this purpose.
In this paper, we first present a heuristic analysis and overview of the ES method, with particular attention to modifications which we find to be advantageous to the algorithm, and to the considerations involved in frequency optimisation. We then present the successful application of an ES controller to the open-loop control system of Oxlade et al. (2015): a highly turbulent bluff-body wake forced by a harmonic pulsed jet, variable in both amplitude and frequency. We demonstrate that in the presence of broadband noise, we are able to optimise both the forcing amplitude and frequency of the open-loop system, either individually or simultaneously in a stepwise manner. We finally demonstrate the ability of the system to adapt in real time to changing conditions, in this case, varying free-stream velocity.
The paper is structured as follows. An overview of the ES method and a summary of the key parameters are given in Sect. 2.1 before the particular considerations involved in frequency optimisation are discussed in Sect. 2.2. A more detailed analysis of the effect of these parameters is then given in Sect. 2.3 and Sect. 2.4. The experimental implementation is then described in Sect. 3 along with the results. Conclusions are finally given in Sect. 4.

Analysis
Before presenting our experimental results, we first give a brief overview of the standard ES algorithm followed by particular issues and suggested modifications. Figure 1a shows schematically two scenarios for the steadystate mapping between a reference parameter r and an output N(r): the upper displaying a desired location with finite gradient and the lower a desired location with zero gradient. Here, we look predominantly at the ES method (for zero gradient) although the same principles and stability guarantees can be applied to both (Ariyur and Krstić 2003). While many variants of ES controllers exist (Calli et al. 2012), the perturbation method is the most commonly applied and is straightforward to analyse and understand. A block diagram for the method is given in Fig. 1b.

Principle of operation
The perturbation method involves adding a perturbation or "dither" signal d(t) to the plant input, the form of which may vary (Nešić et al. 2006;Nešić 2009), but the most commonly used being a sinusoid. The ES system involves an estimation phase, the purpose of which is to produce a signal e, proportional to the local gradient of the mapping, while attenuating the disturbances arising from both the unmodelled dynamics of the plant and the dither signal. The control must act on this estimation signal to produce a change to the reference with the goal of driving the gradient to zero. We will consider first the case of a SISO (single input single output) plant with a nonlinear mapping N(r)(N : R → R); internally generated noise w, and frequency response defined by G(s), where s is the Laplace transform variable. It is important to note here that neither N(r) or G(s) need to be known in order to implement the controller and that G(s) may in general be nonlinear, i.e. amplitude or time dependent. The only assumption is that an operating point can be decided based upon the gradient The equations governing an ES controller applied to this system are given below, while the influence of each of the parameters is summarised in Table 1. Each of the terms featured below may be referred to on the block diagram of Fig. 1b.
The reference r is composed as the sum of some initial value r 0 , the dynamically varying control correction r and the dither signal: As given in Eq. 2 below, the output y from the plant then consists of: a slowly varying part, ȳ; a sinusoidal term proportional to the local gradient of the mapping N ′ (r 0 +r); and noise w G arising from the unmodelled dynamics within (1) r = r 0 +r + a sin(ω |rmd t) .
(a) (b) Fig. 1 a Examples of common nonlinear mappings and (b) the standard extremum seeking control algorithm applied to a SISO system  G(iω), the signal at output will be sinusoidal with a phase lag φ G(iω) with respect to the dither signal and will be attenuated (or amplified) according to a gain |G(iω)| . Such an effect may arise due to the finite response time of the system or due to effects such as hysteresis (Benard et al. 2009).
The high-pass filter F H (iω) then removes the slowly varying term, induces an additional phase lag φ F H and slightly modifies the noise, This output is then multiplied by a unity magnitude dither signal in order to achieve a demodulated signal γ. This operation results in three components: a DC term proportional to the local gradient; a harmonic of the dither signal and modulated noise.
A low-pass filter is then applied to remove the harmonic terms and leave a value proportional to the gradient of the mapping. This will be referred to as the error signal, e, where As will be discussed below, for the case of frequency optimisation, it may be advantageous to choose the dither signal to be a square wave. In this case, the above analysis shows (see "Appendix") that e may instead be written as the Fourier series: Finally, the controller evaluates the correction based upon the value of e with the objective of moving the reference in the direction of the gradient of the mapping. Specifically, K will increase r if e is positive and vice versa: (7) r = K(s)e.
Typically K is chosen to be an integral controller (k / s) in order to drive the system to the desired reference gradient.

Frequency optimisation
As discussed in Sect. 1, for many open-loop flow control systems, the input to the plant may be a harmonic signal, u = c sin (2πft). For example, the input may be a sinusoidally pulsed jet, the effect of which may vary with frequency, as illustrated schematically in Fig. 2a. In many such cases, we may therefore wish to adapt the frequency of the harmonic signal with the ES controller, i.e. f will be the reference parameter r.
Following the standard method described in Sect. 2.1 and applying a sinusoidal dither signal to perturb the frequency, the input to the plant would be, To analyse the effectiveness of this perturbation strategy, we wish to determine the instantaneous frequency content of this signal, since, to determine the gradient of the mapping in For a signal written in terms of a sin function the instantaneous frequency is the derivative of the argument of the sin. That is, for u = c sin(φ(t)), we need to find f (t) = 1 2πφ (t). For the system in Eq. 8 this gives, Hence, the instantaneous frequency is given by the slowly varying base frequency ( f 0 +f ), a sinusoidal perturbation and an additional time increasing term, ω d a cos (ω d t)t. This final term will cause the input to have a far more broadband spectrum, as indicated schematically in Fig. 2b, rather than one that probes the system within the required range. Alternatively, this final term may be thought of as an additional noise component of the reference parameter r. This noise grows with time and is unbounded as t → ∞. This unbounded frequency input will cause erratic behaviour of the system and acts to make the system unstable, in the sense that the reference parameter grows without limit as time increases. Two methods to deal with this problem are now suggested.

Sinusoidal phase deviation
The first method can be found by explicitly calculating the phase argument φ(t) that gives the required variation in frequency, f(t). For a sinusoidal perturbation of the frequency of u = c sin(φ(t)) we require Hence, Here, the approximation results from the assumption that compared with the sinusoid, f changes only very slowly so can be considered constant within the integral. Hence for a sufficiently slow adaptation of f , a sinusoidal dither can be implemented when the reference parameter is a frequency.

Stepwise changes in r
A limitation of the above method comes from the fact that changes in f must be sufficiently slow to avoid introducing undesired frequencies into the input. An alternative solution to the problem is to instead change the forcing frequency in a stepwise manner, akin to the square wave dither signal applied by Nešić et al. (2006). In this case, we can choose the value of f to be updated only when step changes in the dither signal occur, thus avoiding the above problem. The input is now given by (11) Here, we use the subscript T 2k to denote the value at the beginning of each half cycle of the square wave. We can now evaluate the instantaneous frequency at all times other than when the step changes occur as The instantaneous frequency now simply follows stepwise changes with time and will, over the dither period, have the frequency content shown in Fig. 2b. Furthermore, it can be shown (see "Appendix") that this stepwise dither signal will operate in the same way as the sinusoidal dither, giving the result shown in Eq. 6. Therefore, a frequency optimisation scheme can be set up in this way and may then be expected to be subject to the same considerations as the standard ES scheme.

Stability
While proofs of stability for a sinusoidal dither are given in Krstić and Wang (2000) and Krstić (2000) among others, we aim to give a brief heuristic overview here. Instability of the ES system can at first be separated into instabilities of the plant itself and instabilities of the controller. The first situation, instability of the plant, is the most obvious but need not be considered. The purpose of an ES controller is not to stabilise unstable dynamics but to change the average, long term output of the plant. Hence, it should only be used in situations in which the plant itself is open-loop stable. Stability of the controller can now be considered on the basis of the analysis above. The key instance in which the controller can induce instability is if the error signal, e, is of the wrong sign. In general this will be due to phase lags between the dither input and the multiplicative stage. This is clear from the expression for e in Eq. 5. If the total phase lag φ G + φ F H is greater than π/2, the cosine term will change sign, resulting in the control action being in the wrong sense. That is, r will be adapted in the direction away from the peak. Neglecting φ F H this could, for example, happen for a second order plant perturbed above its resonant frequency. An instability of this kind will result in a diverging output from the plant in the opposite direction of the extremum.
Instability may also arise under conditions for which the above analysis is invalid. In particular, it is assumed that the timescales of the dither signal and the adaptation are sufficiently separate. This scale separation is dependent, other things being equal, on the relative values of the control gains K(s) and ω d . A similar observation is noted in the analysis of Krstić (2000) and Tan et al. (2010).
The implication of this limitation will be discussed in Sect. 2.4.

Performance
The performance of the ES system may be considered in terms of three key features, the gradient estimation, disturbance attenuation and control. The performance of estimation and control may be considered separately, as shown schematically in Fig. 1b, although disturbance attenuation is affected by both parts of the system. A summary of the implications of the analysis below is given in Table 1.

Estimation
As shown in Eq. 5, e is composed of the true estimation value, proportional to N ′ (r), and additive noise. For brevity, we will therefore decompose the error signal as e = e t + w e . Estimation quality can now be considered in terms of a signal-to-noise ratio (SNR) for e, defined as, It is shown from Eq. 5 that e t depends upon the amplitude of the dither signal, the local gradient of the mapping and the frequency response of the plant at the perturbation frequency. Of these, only the dither amplitude a can be directly chosen. The local gradient is a property of the plant, is generally unknown, and will vary during operation of the ES system. The plant frequency response may also be uncertain, but its influence may be changed by the choice of ω d . In general, we wish to maximise SNR, hence we wish to maximise the combined term |G| cos(φ G + φ F H ) . It is important to note that all the constituents of this term are frequency dependent, since in general, the magnitude of the frequency response of a physical system decays at higher frequencies, while the phase angle decreases. This consideration therefore provides a limit on the perturbation frequency. The effect of the phase lag may be improved upon by applying a similar phase lag to the demodulating signal, that is, to multiply by sin(ω d t − ϕ) (Krstić 2000). If ϕ can be chosen to be very similar to φ G then the SNR may be improved; however, in general φ G may be unknown and limitations on ω d will still remain.

Disturbance attenuation
The other approach to maximising the SNR is to minimise the noise, w e . A significant part of this noise comes from the dither signal and its harmonics: demodulation leads to the sinusoidal term in Eq. 4, while nonlinearities within the plant may lead to higher harmonics (e.g. 2ω d , 3ω d ,...) of the dither signal in the output y. These harmonic terms can be especially damaging to the control performance because, being at relatively low frequencies of O(ω d ), they are attenuated very little by an integral controller. Successful elimination of these terms therefore depends upon the properties of F L (iω). A particularly good choice for this filter is a moving average filter with period T = 2π/ω d , since the frequency response of such a filter has zero gain at multiples of ω. This is shown in Fig. 3. In general, the majority of the remainder of the noise will originate from the internal dynamics of the plant, and will have spectral characteristics that depend on the plant properties. For example, these spectral properties may be similar to those of the plant transfer function G(jω). In general, these noise components will be attenuated by the moving average filter as well as by the integral controller.
While the low-pass filter is intended to remove noise after the demodulation, the purpose of the high-pass filter is to remove the DC component ȳ from y, thereby isolating the harmonic perturbation. Exclusion of the high-pass filter simply leads to an additional harmonic term in γ (see Eq. 4) of frequecy ω d . However, as discussed, such a harmonic may be effectively removed by the low-pass filter. The high-pass filter may therefore be redundant. In summary, we may expect to be able to remove the high-pass filter, thereby avoiding the associated phase lag. Although this lag may be catered for by including an additional lag in the demodulating signal, it may still impose speed restrictions on the adaptation of the whole algorithm.
Control Historically, ES controllers have often used a pure integral controller (Ariyur and Krstić 2003), K(s) = k/s. This guarantees zero steady-state error and has the added benefit of acting as a low-pass filter to the many disturbances passing around the loop. However, if the gradient estimation is of sufficient quality, then K(s) can Fig. 3 Frequency response of the moving average filter of period T = 2π/ω d be designed to be more aggressive, using a PID or more general frequency domain controller. For example, Krstić (2000) proposes that K(s) include a dynamic compensator, the properties of which are decided based upon G(s) and the properties of the expected changes in the operating point.

Experimental Implementations
Based upon the analysis above, we apply an ES controller to the open-loop control system of Oxlade et al. (2015). The open-loop system here has a steady-state mapping between the parameters of sinusoidal pulsed jet forcing and the base pressure (corresponding to the drag) of an axisymmetric bluff body. For this system, the mapping is two-dimensional; mean base pressure is a function of both forcing amplitude and forcing frequency. The experimental setup and system properties will now be described in more detail.

Experimental setup
The principal features of the experimental setup are the same as those described in Oxlade (2013) and Oxlade et al. (2015) consisting of an axisymmetric, bullet-shaped bluff body fitted with a pulsed jet actuator. A schematic of the experiment is shown in Fig. 4. The model is instrumented with 8 Endevco 8501C-1 pressure transducers, and 64 static taps, distributed on a polar grid as shown in Fig. 4. The pulsed jet actuation is applied by means of a speaker located within a cavity inside the rear of the model. Oscillation of the speaker diaphragm forces air through an annular slit on the back face of the body, the resulting zero-netmass-flux jet emanating in a stream-wise direction.
The model was tested in a closed-loop tunnel operating with a free-stream velocity of 15 ms −1 , giving a Reynolds number based on body diameter of Re D = 1.88 × 10 5 , and a fully turbulent wake.The Endevco transducers were used during real-time implementations of the ES algorithm, while the 64 static tappings were used to calculate the mean base pressure in post-processing. The free-stream velocity was measured via a Pitot-static tube and the tunnel speed PID controlled to achieve constant velocity.
Data acquisition and implementation of the controller were both performed on a National Instruments PXIe-1078 chassis with a PXIe-6358 data acquisition card, and running the real-time operating system. The control loop and data acquisition were performed at 2 kHz, while the output to the speaker was sent to the digital to analogue converter of the acquisition card at 90 kHz. While the frequency of the forcing could be controlled directly from the actuation signal, forcing amplitude was quantified in terms of the root mean square (rms) pressure ( p 2 c ) within the cavity, measured via an additional pressure transducer (p c ). We use the rms cavity pressure as this scales with the peak jet velocity (Oxlade 2013) and provides the most practical measurement for use in real-time applications. A calibration between p 2 c and the jet velocity was not performed here as this is specific to the free-stream velocity and, moreover, is not required for the control algorithm. Control of the rms pressure was implemented via a PID loop in order to cater for the frequency response of the pulsed jet system and for any changes in operating conditions.

System properties
For an ES controller, the most important feature of the plant is the steady-state mapping N(r). For the system in this study, the mapping is two-dimensional (N : R 2 → R ), giving spatially averaged base pressure as a function of forcing amplitude and frequency as shown in Fig. 5. Here, the temporally and spatially averaged base pressure change is defined as where A is the area over the base, Ĉ p is the pressure coefficient without forcing and C p is the pressure coefficient under forced conditions. The mapping is smooth with a single global maximum located at values of amplitude and frequency of 1200 Pa rms and 750 Hz, respectively. Areas for which data are not shown are outside the range of the actuator and therefore cannot be explored in open-loop or reached during operation of the ES controller.
For a bluff body such as that tested here, much of the drag is form drag. A positive C p therefore corresponds to a drag reduction. The mapping displays a drag reduction for large amplitudes and for frequencies f 500 Hz, corresponding to a situation in which the forcing frequency is sufficiently decoupled from those of the shear layer and coherent structures of the wake. The forcing jet has zeronet-mass-flux, but is not referred to here as a synthetic jet as the mechanism of operation does not rely on self induced streaming. Instead the forcing generates a discrete train of vortices bounded by strong shear layers that provide a sheltering effect on the wake, reducing entrainment. We refer the reader to Oxlade et al. (2015) for further details on the drag reduction mechanism and actuator details.
While the data shown here is for a fixed free-stream velocity U ∞ = 15 ms −1 , one might expect the mapping to evolve in a predictable way as conditions change. Dimensional analysis and knowledge of the drag reduction mechanism tell us that C p is a function of jet momen- Here, u j and A j are, respectively, the jet velocity and area, f is the forcing frequency and θ is the boundary layer momentum thickness at separation. Provided that the Reynolds dependence of the dimensionless mapping is small, it may be expected that the optimal jet velocity and frequency will scale linearly with U ∞ .

Extremum seeking algorithm
The implemented ES controller was operated to adjust the amplitude and frequency of the harmonically forced jet described above, either of these parameters thereby acting as the reference r described in Sect. 2.1. The output of the system (y) is taken to be C p , approximated by the average of the eight pressure transducers and evaluated independently for each time sample.
The implemented ES algorithm is shown in the block diagram of Fig. 6. Following our analysis in Sect. 2.4, the implemented system had the following features not commonly used in the literature:  Oxlade et al. (2015). Spatially averaged base pressure is given as a function of forcing amplitude and frequency, where the forcing amplitude is defined in terms of the rms of the pressure (p c ) inside the cavity. Negative contours are shown by the dashed lines Fig. 6 Modified ES algorithm implemented experimentally. The output y is the pressure instantaneously averaged over the base, while the reference r is either the forcing amplitude or frequency 1. Square-wave dither signal for frequency As discussed in Sect. 2.2, a sinusoidal dither signal leads to complications when the reference parameter is itself a frequency. We therefore choose to implement the method described above, varying the forcing frequency in a stepwise manner, and only updating r at the crossing points of the square wave.

Exclusion of the high-pass filter
As discussed in Sect. 2.4, the main purpose of the highpass filter is to remove the DC component of the plant output, thereby removing a sinusoidal term generated by the demodulation. Provided that this effect can be successfully achieved using the low-pass filter, the high-pass filter is unnecessary. Furthermore, the highpass filter introduces an additional phase lag which may impose speed restrictions.

Use of a moving average for the low-pass filter
The key components that need be removed by the lowpass filter are the perturbation frequency and its harmonics. A moving average filter with period equal to the perturbation period will have zero gain at all harmonics of the inverse of this period (Smith 1997), and is therefore the optimal filter to use for this purpose (see Fig. 3).

Use of a proportional, integral controller
To improve the speed of adaptation, a proportional term is used in addition to the usual integral term. This is in agreement with the recommendations of Krstić (2000).

Parametric investigation
Given the overall structure of the ES system described above and shown in Fig. 6, it remains only to choose the parameters of the system. If the ratio of proportional to integral gains is fixed, then only three parameters remain to be chosen: the dither amplitude a, dither frequency ω d and control gain k. These are defined such that The ratio of the proportional and integral gains shown here was based on the recommendations of Krstić (2000) and some initial trial and error. The results of the parametric investigations for seeking in amplitude are shown in Table 2. A preliminary investigation was first performed to establish approximately optimal values for these variables of a = 100 Pa, ω d /2π = 0.5 Hz and k = 0.8. Subsequently, these initial values formed a baseline about which the parametric investigation could be performed; i.e. for investigation of a, ω d /2π is kept at 0.5 Hz and the control gain k kept at 0.8. The table displays the half rise time, the steady-state variance and the stability of the scheme under each configuration. The half rise time is defined as the time taken for the reference to reach half of the average steady-state value, thereby giving a measure of the convergence speed. The steady-state variance is the variance of the reference value after it has reached a steady-state, including both the dither signal and the control input r. The system is determined to be unstable if the reference parameter increases or decreases beyond the range of the actuator and well beyond the anticipated optimal condition.
In agreement with the analysis of Sect. 2.1, the dither amplitude is seen to have two key influences. An increase in a is seen to increase the speed of the algorithmthe rise time reducing by a factor of 10 with a threefold increase-but is also naturally seen to increase the perturbations seen once a steady state is reached. Within the range tested here, the system always remained stable. By contrast, the dither frequency is seen to have little coherent effect on the convergence time, although at very low frequencies there seems to be more perturbation about the steady state. At high frequencies, the system is found to become unstable, as anticipated: the reference increasing beyond the optimal operating point. This may be a result of the phase lag through the system becoming too large. Finally, an increase in the control gain is seen to increase the convergence speed, as anticipated, but is also seen to lead to greater noise in steady state due to overreaction of the controller.
Choice of optimal parameters require a choice of tradeoff between convergence speed and steady-state perturbations. It was felt that the values of a = 100 Pa, ω d /2π =  )  -46  63  78  101  99  89  77  81  --73  78 108 0.5 Hz and k = 0.8 provided a good compromise and were used for all subsequent tests. Similar results for frequency optimisation lead to identical choices for ω d and k and a dither amplitude a = 30 Hz.

Single variable implementations
The ES controller was first implemented in one dimension to optimise either only the amplitude of the forcing at a fixed forcing frequency or to optimise the frequency at a fixed amplitude. Initial conditions for the amplitude optimisation were 400 Pa and 650 Hz, while those for the frequency optimisation were 250 Pa and 200 Hz, chosen in order to provide sufficient "space" in which the controller could seek. Examples of these implementations are shown in Figs. 7 and 8.
For the case of amplitude optimisation, the response of the output of the system under control action is shown in Fig. 7a. The trajectory of the system across the static mapping is shown inset. Two key observations can be made from these data. First, by looking at the filtered response (black line), the output can be seen to have adapted to an optimal condition after around 30 s, beyond which only fairly small fluctuations occur. Second, the raw data (blue) displays the level of noise resulting from the broadband turbulent fluctuations of the flow, all of which pass into the output measurement C p . These fluctuations can be seen to have been greatly attenuated within the error signal e, shown in Fig. 7b, although some fluctuations do still remain. Finally, the adapted reference input can be seen as the black line in Fig. 7b. It is evident that the reference has been largely adapted after 30-40 s and that only small fluctuations remain after this time. The filters and controller  are therefore effective in deciphering the effect of the dither within the noisy output signal y, and passing only minimal disturbances into the adapted reference r.
Similar results can be seen for the case of frequency optimisation in Fig. 8. Adaptation is seen to be similarly effective in spite of the same noisy output measurement. However, the key difference relative to amplitude optimisation is that the frequency is seen to continue to increase with time, because the gradient of the mapping remains slightly positive. For this particular operating point, a slope (rather than extremum)-seeking controller may therefore be more appropriate; however, for the global optimum, this is not the case. Also shown is that changes in the frequency are made in the stepwise manner described in Sect. 2.2.
The adaptation time for both controllers is seen to compare well with previous flow control implementations. In particular, the adaptation time is found to be similar to, but slightly greater than, that of the system of Henning et al. (2008) and Pastoor et al. (2008). In general, for a given noise level from the plant (w of Fig. 1), there exists a tradeoff between the adaptation time and the level of variation once a steady state is reached, as shown in the parametric study of Sect. 3.4. For a turbulent wake, we may expect the level of noise to be partly determined by the Reynolds number as with increasing Re there is an increased range of scales in the flow. If we wish to keep the steady-state perturbations below a certain level, the adaptation time would be required to increase with Re. The slight increase in adaptation time may therefore be at least partially attributed to the factor 10 higher Reynolds number in this experiment.

Dual variable implementation
While the implementation of the algorithm in each of amplitude and frequency individually allowed optimisation of the algorithm parameters, a fully working system would be operated in both dimensions. In principle, the ES controller can be operated to adapt multiple references simultaneously, provided that the perturbation frequencies are not equal; a result of the orthogonality of sinusoids. However, the challenge arises in effectively filtering the variable γ (see Fig. 6). The moving average filter discussed in Sect. 2.4 is very effective at removing one frequency and its harmonics but not optimal for other frequencies.
Each moving average is therefore unable to effectively filter the dither signal applied to the other variable, since the frequency of the two dither signals must necessarily be different.
While methods do exist to avoid the filtering issue (see e.g. Gelbert et al. 2012), in this study, we chose to implement the algorithm in a stepwise manner, adapting each variable in turn for a fixed period of time. The results of this are shown in Fig. 9. The system converges to a steady state after about 100 s, so is clearly somewhat slower than the one-dimensional examples seen above. This may be partially a result of the stepwise nature of the control, as the controller has to effectively restart from its current position every 20 s. However, the total convergence time still compares favourably with previous studies (Pastoor et al. 2008).
The trajectory of the adaptation can be seen in Fig. 9b, from which it is clear that the system converges to the optimal condition found in open-loop. Each stage of the route is seen to arrive at approximately optimal conditions based upon the current value of the fixed parameter (amplitude or frequency). It is possible that by reducing the time between switches from the 20 s implemented here, the net speed of the algorithm may be improved. However, the best time will likely depend on the local properties of the mapping so a generally optimal value may be difficult to find. It is worth noting that under these testing conditions, the controller brings the system close to the limits of the actuator. While this is not an issue to the control scheme here, it is clear that for any practical system, care must be taken to ensure that the controller does not try to reach an operating point that would cause damage to the actuator.

Adaptation to changing conditions
One of the key purposes of using closed-loop control over open-loop control is that the system will automatically adjust to changing operating conditions. For this ES control system, the key operating condition is the free-stream velocity U ∞ , since this affects the optimal amplitude and frequency as discussed in Sect. 3.2. Specifically, with increasing U ∞ , the required forcing amplitude and frequency would be expected to increase accordingly. The system was therefore tested with a slow sinusoidal variation in U ∞ , optimising in forcing amplitude only. The results of this are shown in Fig. 10.
The results indicate that the controller is able to adapt to the changing conditions well and with minimised disturbances. Over the first 50 s, the system converges to the optimal value of p 2 c ≈ 1000 Pa, corresponding to U ∞ = 15 ms −1 . With varying velocity between 10-20 ms −1 , the amplitude is seen to also vary between approximately 800-1600 Pa with a small lag. Adaptation is therefore larger in the positive direction than in the negative direction, possible indicative of nonlinearity in the mapping between p 2 c and u j . Regardless, the controller is here demonstrated to be able to adapt in both a positive and a negative direction.
While it is not possible to confirm that an optimal condition is maintained at all times, a strong adaptation is evident. While truly optimal adaptation must take place in frequency as well as amplitude, the results of shown in Fig. 9 indicate that the controller would be able to achieve this, provided that variations are sufficiently slow.

Concluding remarks
We have provided a heuristic analysis of the ES algorithm with a particular focus on flow control applications. Our analysis has demonstrated that the often used high-pass filter may be removed, and that the low-pass filter may be replaced with a moving average in order to optimally filter the perturbations resulting from demodulation. It has also been demonstrated that for the case where the reference parameter is the frequency of a harmonic input, implementation of a sinusoidal dither signal is problematic. In this case, a stepwise dither signal can provide a suitable alternative, provided that the control adaptation r is only updated when the step changes in the dither signal occur. This therefore provides a suitable ES method for a number of open-loop flow control applications involving a harmonic input.
Based on the analysis of the algorithm, a modified ES controller was implemented experimentally in order to demonstrate its efficacy. The controller was applied to the open-loop control system of Oxlade et al. (2015), providing a closed-loop extension to the pre-existing system. Despite a noisy output signal perturbed by turbulent fluctuations, the controller is able to find optimum forcing configurations within around 30-40 s when operated on a single forcing variable and around 100 s when operated on two forcing variables. Furthermore the system is able to effectively filter the measured noise rather than passing it to the adaptation signal, and is therefore able to maintain a fairly steady output once convergence has taken place. The controller is also able to adapt to changing conditions, in this case free-stream velocity, in real-time. While in this investigation, the controller was implemented to seek maximum drag reduction, the system could just as easily be used to find an optimal condition for energy efficiency, as in Beaudoin et al. (2006).
In conclusion, we have introduced a modified ES algorithm and demonstrated its efficacy in adjusting the openloop forcing of a turbulent bluff-body wake. We hope that the modified system proposed in this paper will be of use for other flow control applications, particularly those involving a harmonic input signal.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix: square-wave dither signals
In Sect. 2.2, it was proposed that a square-wave dither signal could be used to overcome the issue of varying the frequency of an input signal. To check that a square-wave perturbation and associated square-wave demodulation will work, we can write the signal in terms of its Fourier series. For a square wave s(t) of period T = 2π/ω d ,
The items within this summation can be split into two cases after applying the trigonometric identity For m = n we have a sum over terms of the form given in Eq. 4 which include a DC component, while for m � = n we have purely sinusoidal terms. Given a sufficiently effective low-pass filter, only the DC terms will pass into e: This square-wave ES system can therefore be expected to work in a very similar manner to the original system with a sinusoidal dither, and is suitable for use in optimising the frequency of a harmonic signal.