2.1. Aerofoil model and plasma actuator

Experiments were carried out in the open-circuit low-speed wind tunnel at the National Key Laboratory of Aerospace Power System and Plasma Technology. As shown in figure 1, the tunnel has a 2-metre long square test section, with a cross-sectional area of 0.2 $\times$ 0.2 m. The test section is made of an acrylic plate (thickness: 20 mm) for optical imaging. The maximum free-stream velocity in the test section can reach 50 ms⁻¹ and the typical turbulence level is below 0.2 %. A classical symmetrical aerofoil (NACA-0015) three-dimensionally printed in resin is installed vertically in the test section, and the angle of attack ( $\alpha$ ) could be adjusted by a rotary platform. The chord length (c) and span (h) of the aerofoil are 0.15 m and 0.18 m, respectively. In this study the blockage ratio of the aerofoil model calculated from $\varepsilon =ct/A_{\mathit{t}}$ reaches $8.44\,\%$ , where $t$ is the maximum thickness of the aerofoil and $A_{\mathit{t}}$ is the cross-sectional area of the test section. To avoid overestimation of the aerodynamic coefficients, the blockage correction method proposed by Maskell (Reference Maskell1972) is employed (see Appendix B for details).

Figure 1.

Sketch of the aerofoil model in the test section (not to scale): (a) top view, (b) side view.

Figure 1. Long description

Panel A: Top view of the aerofoil model. The diagram shows a freestream flowing from left to right. The aerofoil is mounted on a pivot and can rotate around an axis labeled as alpha. A hotwire sensor is positioned downstream of the aerofoil to measure flow characteristics. A step motor is connected to the side wall to control the position of the aerofoil. The field of view is indicated by a rectangular area. Panel B: Side view of the aerofoil model. The aerofoil is positioned vertically within a test section. A plasma actuator is mounted on the aerofoil, and a laser beam is directed towards it. A hotwire monitor station is placed to the right of the aerofoil. The setup includes a rotary platform and a force sensor at the base to measure forces exerted on the aerofoil.

The bottom of the aerofoil model is attached directly to a two-component force sensor (Autoda 9102, range: 0–10 $N$ , measurement uncertainty: 0.1 %). To get rid of the influence of the wall boundary layer, the two ends of the aerofoil are 10 mm away from the floor and ceiling of the test section. In this study the free-stream velocity ( $U_{\infty }$ ) is fixed at 15 ms⁻¹. Based on this velocity scale and the aerofoil chord length, the Reynolds number is determined to be ${\textit{Re}}_c=1.5\times 10^5$ . A coordinate system is established at the middle point of the aerofoil leading edge, with the x, y and zaxes along the streamwise, wall-normal and spanwise directions, respectively. The hotwire sensor is placed in the symmetry plane ( $z=0$ mm), with its streamwise coordinate fixed at $x=1.01c$ and its spanwise position adjusted by a stepper motor.

Figure 2.

(a) Structure of the DBD plasma actuator. (b) Representative discharge voltage and current waveforms.

To mitigate the flow separation at high angles of attack, a DBD plasma actuator is flush mounted on the suction surface of the aerofoil, at approximately 3 % chord length away from the leading edge (see figure 1 a). As shown in figure 2(a), the plasma actuator consists of a high-voltage copper electrode (thickness: 0.05 mm, width: 3 mm), a Kapton dielectric layer (thickness: 0.18 mm) and a ground copper electrode (thickness: 0.05 mm, width: 5 mm). The length of the actuator is 180 mm, covering the entire aerofoil span. High-voltage wires are welded to the copper electrodes for subsequent electrical connection with the power supply. A high-voltage sinusoidal waveform generator (GBS Elektronik GmbH, Minipulse 2.2) is used to drive the plasma actuator, and the output voltage amplitude scales proportionally with the input voltage level ( $V_{d}$ , range: 12–40 V). Two digital signals are required to fully define the output waveform. The first signal is a square-wave clock signal that controls the base frequency ( $f_b$ ) of the sinusoidal waveform, and the second signal is an enable signal that determines the on/off status of the plasma actuator. In open-loop periodical control the enable signal is set as a periodical pulse signal, with adjustable duty cycle ( $D_c$ ) and modulation frequency ( $f_m$ ). For the purpose of present study, $f_b$ and $V_{d}$ were fixed as 10 kHz and 26 V, respectively.

A high-voltage probe (Tektronix P6015A) and a current monitor (Tektronix P6022) are attached to the anode and the ground lead to measure the discharge voltage ( $U_{d}$ ) and the discharge current ( $I_{d}$ ), respectively. Figure 2(b) shows the representative discharge waveform of the plasma actuator in steady working mode (i.e. $D_c=100$ %). At a peak-to-peak discharge voltage of 9 kV, the maximum discharge current amplitude reaches about 150 mA, occurring at the zero crossing point of the rising stage of the voltage waveform. By integrating the voltage and current within a discharge period (denoted as $T_d$ ), the average power consumption of the actuator ( $P_{\mathit{act}}$ ) is estimated to be 4.7 W in steady working mode, as shown in (2.1):

(2.1) \begin{equation} P_{\mathit{act}}=\frac {1}{T_d}\int _{0}^{T_d}U_dI_d\boldsymbol{\cdot }\text{d}t. \end{equation}

In studies related to drag reduction, a dimensionless metric called control efficiency is typically used to measure the true performance of control, which can be defined as the ratio of the energy saved by drag reduction to the energy consumed by plasma actuation. For the present study, the control objective is to maximise the lift coefficient of the aerofoil model, and a proper definition of the control efficiency is not straightforward. Following the principle of dimension matching, the following definition of the control efficiency is adopted:

(2.2) \begin{equation} \eta =\Delta F_l \boldsymbol{\cdot }U_{\infty }/P_{\mathit{act}}. \end{equation}

Here $\Delta F_l$ denotes the lift increment caused by plasma actuation.

2.2. Control and measurement system

Figure 3 shows the signal flow in the control and measurement system. Three core components constitute the closed-loop control system, including a real-time controller (National Instruments, CRIO-9038), a hotwire anemometer (Hanghua Technology, CTA-04) and a DBD plasma actuator. In the chassis of the controller, two data acquisition modules are installed: an analogue input module (NI-9223) for high-speed hotwire signal acquisition and a digital output module (NI-9264) for rapid plasma actuation adjustment. Both of them are hooked directly to a high-performance FPGA chip (Kintex-7 160T) embedded in the controller, with an on-board clock frequency reaching as high as 40 MHz, which enables real-time closed-loop control at the microsecond time scale. In Rabault & Kuhnle (Reference Rabault, Kuchta, Jensen, Réglade and Cerardi2019) the influence of control frequency ( $f_c$ ) on the long-term reward of DRL was investigated via numerical simulation. They found that when $f_c$ was much lower than the characteristic frequency of the flow ( $f_0$ ), the DRL agent failed to learn effective control strategies because the update of control action is too slow to catch up with the rapidly evolving fluid dynamics. On the contrary, if $f_c$ was too high (say, two orders higher than $f_0$ ), the learning also failed because the agent could not feel any impacts of its action in the near future. Due to the above concern, $f_c$ is selected as 1000 Hz in the present study, one order higher than the characteristic frequency of the aerofoil flow: $f_0=U_\infty /c=100 \,\mathrm{Hz}$ . For training purpose, the real-time controller communicates with a host computer via a network switch. The model numbers of the central processing unit (CPU) and graphics processing unit in the host computer are Intel Core i9-12900E and NVIDIA GeForce RTX 3060, respectively.

Figure 3.

Signal flow in the control and measurement system.

As shown on the right-hand side of figure 3, three types of measurements are conducted, including single-point hotwire measurement, direct force measurement and particle imaging velocimetry (PIV) measurement. In Zong et al. (Reference Zong, Wu, Li, Su and Liang2024) a closed-loop control of the aerofoil trailing-edge flow separation was implemented with Q-learning, and an anti-correlation was found between the instantaneous flow separation area and the velocity measured by a hotwire at the separated shear layer ( $U_h$ ). Inspired by the above study, the hotwire sensor was positioned in the aerofoil wake zone to provide a real-time feedback of the control effect. The y position of the hotwire (denoted as $y_h$ ) can be adjusted conveniently by a step motor. During hotwire calibration, the aerofoil is removed and the hotwire sensor is placed in the mainstream. In total, 10 calibration points are arranged within the range of 0–25 ms⁻¹ and, for each point, the output voltage of the hotwire is sampled for 10 s at a sampling rate of 10 kHz. The relationship between the mean hotwire voltage and the wind speed is fitted by a four-order polynomial, with a typical fitting error below 0.5 %. To compensate for air temperature variations in the wind tunnel, the temperature-correction calibration method developed by Hultmark & Smits (Reference Hultmark and Smits2010) is implemented.

In this study, high-voltage discharge inevitably introduces electromagnetic interference to the control system, and several measures are taken to minimise its adverse impact. First, the high-voltage power supply and the data acquisition devices are fed by two independent electrical systems, physically separated by an isolation transformer, to avoid the electrical noise introduced from common ground. Second, the signal transmission cables connecting the hotwire anemometer, the data acquisition device and the real-time controller are wrapped with two layers of fine copper mesh to form a Faraday cage, which shields the noise introduced by spatial electromagnetic waves. Third, after the raw hotwire voltage is sampled, a median filter (length: 3) is used to remove the abnormal spikes in the measured signal.

The aerodynamic forces of the aerofoil model are diagnosed by a two-component force sensor attached at the bottom (see figure 1). Since the force sensor rotates simultaneously with the aerofoil, the measured two force components (denoted as $F_x$ and $F_y$ ) are always along with the chord direction and wall-normal direction, respectively. The lift coefficient ( $C_l$ ) and drag coefficient ( $C_d$ ) are related with $F_c$ and $F_n$ as

(2.3) \begin{equation} \left \lbrace\!\! \begin{array}{ll} C_l=\displaystyle \frac {F_y\cos \alpha -F_x\sin \alpha }{0.5\rho _0 U_{\infty }^{2}\boldsymbol{\cdot }c\boldsymbol{\cdot }h}, \\[12pt] C_d=\displaystyle \frac {F_y\sin \alpha +F_x\cos \alpha }{0.5\rho _0 U_{\infty }^{2}\boldsymbol{\cdot }c\boldsymbol{\cdot }h} ,\\ \end{array} \right . \end{equation}

where $\rho _0$ is the air density.

The flow field above the aerofoil suction surface in the $xy$ plane is measured by a high-speed planar PIV system, which consists of a high-speed camera (Phantom VEO-310L, resolution: 1280 $\times$ 800 $\mathrm{pixels}^2$ ), a high-speed laser (Grace Laser, TABOR-D30, pulse energy: 30 mJ) and a programmable timing unit (LaVision, PTU-X). To avoid interference between the hotwire sensor and the laser sheet (thickness: 1 mm), the PIV measurement plane was positioned at a spanwise station of $z = 10\,\mathrm{mm}$ , slightly offset from the hotwire measurement plane ( $z = 0\,\mathrm{mm}$ ). The camera is equipped with an objective of 50 mm (Nikon, Micro Nikkor) and imaging a field of view of 165 mm $\times$ 90 mm (1.1 $c$ $\times$ 0.6 $c$ ). The raw particle images were recorded and processed with Davis 10.2. During post-processing, the final cross-correlation iteration employed an interrogation window size of 32 $\times$ 32 $\mathrm{pixels}^2$ and an overlapping ratio of 75 %, yielding a final spatial resolution of 2.24 mm. Table 1 lists all the cases selected for PIV measurement. In the open-loop control mode, five cases with increasing modulation frequencies of $F_m= 10,\,50,\,100,\,150,\,200,\,400\,\text{Hz}$ are tested. To obtain statistically converged velocity components, an odd sample rate of 1001 Hz is used to avoid having any common divisor with the modulation frequency. For DRL control, both the training stage and evaluation stage are diagnosed. For comparative cases such as deep Q network (DQN) and linear genetic programming (LGP), only the best control law is tested. The evaluation case shares the same acquisition setting as the open-loop control case. In the training process, since it lasts for 50 s, the PIV sampling rate is reduced to 100 Hz. Peak detection error and finite ensemble size are identified as the primary sources of velocity measurement uncertainty. Following the procedure depicted in Zong et al. (Reference Zhu and Krstic2022), the overall measurement error of the mean velocity can be computed to be 1.87 %.

Table 1.

Cases selected for PIV measurements.

Prior to presenting any results, it is necessary to make a convention on the symbol usage. The instantaneous velocity components along the $x$ , $y$ and $z$ axes are denoted as $u$ , $v$ and $w$ , respectively. The time-averaged velocity components along the three axes are represented by $U$ , $V$ and $W$ , respectively. By analogy, $u'$ , $v'$ and $w'$ stand for the three fluctuating velocity components. Ensemble-averaged operation is denoted as angle brackets ’ $\langle \rangle$ ”. Here $\nu$ and $P$ are the kinematic viscosity and the local pressure of the air, respectively; $k$ and $P_k$ stand for the turbulent kinetic energy (TKE) and the turbulence production rate, respectively. Other symbols will be introduced at their first occurrence.

3.1. Force measurements result

Figure 4(a) shows the baseline lift and drag coefficients measured at $U_\infty =$ 15 ms⁻¹. In the range of $\alpha \leqslant 16^\circ$ , the lift and drag coefficients increase linearly with the angle of attack, as prescribed by the thin aerofoil theory (Anderson Reference Anderson2007), the maximum lift coefficient reaches 1.22. At $\alpha = 18$ $^\circ$ , the aerofoil stalls, manifested by a sudden drop of $C_l$ and a sharp rise of $C_d$ . Since our purpose is to reduce the large-scale leading-edge separation, a deep stall angle of $\alpha = 20^\circ$ is selected to perform subsequent flow control experiments. Under open-loop control, the output voltage amplitude is fixed, leaving the modulation frequency ( $F_m$ ) and duty cycle ( $D_c$ ) as two independent variables. In total, five discrete levels of $D_c$ ( $10\,\%,\,30\,\%,\,50\,\%,\,70\,\%,\,100\,\%$ ) and six levels of $F_m$ ( $10,50,100 ,150 ,200 ,400 \,\text{Hz}$ ) are selected to perform the cross-test, resulting in a total of 30 parameter combinations. For each combination, three repeated force measurements are conducted, and the corresponding control effectiveness is quantified by the variation of mean lift coefficient as

(3.1) \begin{equation} \Delta C_l = \frac {1}{3}\sum _{i=1}^3 \left ( C_{l,\textit{AT}}^{(i)} - C_{l,\textit{BL}}^{(i)} \right ), \end{equation}

where $C_{l,\textit{AT}}$ and $C_{l,\text{BL}}$ denote the lift coefficients under plasma actuation and baseline conditions, respectively.

Figure 4.

(a) Lift and drag coefficients under baseline conditions. (b) Plasma-induced lift coefficient increment ( $\Delta C_l$ ) under open-loop control. The blue dots represent the sampling points. The thin dashed black lines are contour lines in a step size of 0.03.

Figure 4. Long description

Panel A: A line graph shows the lift and drag coefficients under baseline conditions. The x-axis represents the angle of attack (alpha) in degrees, ranging from 0 to 20 degrees. The left y-axis represents the lift coefficient (C_l) in blue, ranging from 0 to 1.25. The right y-axis represents the drag coefficient (C_d) in red, ranging from 0 to 0.5. The blue line indicates an increasing trend in the lift coefficient as the angle of attack increases, peaking around 15 degrees. The red line shows a gradual increase in the drag coefficient with the angle of attack. Panel B: A contour plot illustrates the plasma-induced lift coefficient increment (ΔC_l) under open-loop control. The x-axis represents the Strouhal number (S_t), ranging from 0.1 to 4.0. The y-axis represents the duty cycle (D_C) in percent, ranging from 0 to 100 percent. The contour lines, represented by thin dashed black lines, show increments in a step size of 0.03. The color gradient from light to dark red indicates increasing values of ΔC_l, with the highest values concentrated around the middle of the plot. Blue dots represent sampling points.

In figure 4(b) the variations of $\Delta C_l$ at different actuation parameters are plotted as a two-dimensional map. The modulation frequency $F_m$ is normalised by the aerofoil chord length $c$ and free-stream velocity $U_\infty$ , which gives rise to a dimensionless quantity called the Strouhal number, $St = F_m c/U_\infty$ . A noticeable increase of the lift coefficient is observed in the majority of the mapped area, except for the scenario of steady plasma actuation at $D_c=100\,\%$ . The maximum value of $\Delta C_l$ reaches 0.148, located at an electrical setting of $St = 1$ and $D_c = 10$ %. The optimal Strouhal number for the flow separation control is consistent with the review in Greenblatt & Wygnanski (Reference Greenblatt and Wygnanski2000). When the Strouhal number is fixed at 1, the magnitude of $\Delta C_l$ drops monotonically with increasing duty cycle, indicating that the momentum of plasma actuation does not play a crucial role in this flow control scenario.

Figure 5.

Time-averaged velocity fields of the aerofoil under (a) baseline conditions, and (b–f) open-loop control. The duty cycle is fixed at $D_c = 10$ %. In each plot, the in-plane velocity magnitude is shown as a contour. The grey arrows and red lines indicate vertical velocity profiles at representative chordwise locations. The solid black line marks the zero-velocity line ( $U/U_\infty = 0$ ) under actuated conditions. The thin black dashed line indicates the zero-velocity line under baseline conditions.

Figure 5. Long description

The image contains six contour plots labeled (a) through (f), showing the time-averaged velocity fields of an aerofoil under baseline conditions and open-loop control with varying Strouhal numbers. Each plot includes a color gradient representing the in-plane velocity magnitude, with a color bar at the top indicating the velocity ratio ranging from 0 to 1. The grey arrows and red lines indicate vertical velocity profiles at representative chordwise locations. The solid black line marks the zero-velocity line under actuated conditions, while the thin black dashed line indicates the zero-velocity line under baseline conditions. Panel A shows the baseline condition with St = 0. Panel B shows the condition with St = 0.1. Panel C shows the condition with St = 0.5. Panel D shows the condition with St = 1. Panel E shows the condition with St = 2. Panel F shows the condition with St = 4. The x-axis is labeled x/c and the y-axis is labeled y/c, with both axes ranging from 0 to 1 and −0.3 to 0.1, respectively.

3.2. Particle imaging velocimetry diagnosis and discussion

Figure 5 compares the time-averaged velocity fields measured by PIV in baseline and open-loop control conditions. Here $U_{\textit{xy}}$ denotes the in-plane velocity magnitude computed from the Euclidean sum of the streamwise velocity component $U$ and the wall-normal velocity component $V$ :

(3.2) \begin{equation} U_{\textit{xy}}=\sqrt {U^2+V^2}. \end{equation}

Under the baseline scenario of $\alpha = 20^\circ$ , a large-scale leading-edge flow separation is exhibited on the suction side of the aerofoil, with a prominent reverse flow residing in the near-wall region. The contour line of $U_x=0$ ms⁻¹ (i.e. zero-velocity line) separates the forward and reverse flow, defining the time-averaged position of the separated shear layer. After periodical plasma actuation (see figure 5 b–f), although the separation point always remains near the leading edge, a slight reduction of the recirculation area ( $A_{\mathit{b}}$ ) is clearly observable for all cases. Quantitatively, $A_{\mathit{b}}$ is defined as the region bounded by the zero-velocity line, the aerofoil model and the vertical line intersecting the trailing edge. Taking the square of the chord $c$ as the reference area, a dimensionless flow separation area $A_{\mathit{b}}^*$ can be derived as

(3.3) \begin{equation} A_{\mathit{b}}^*=\frac {A_{\mathit{b}}}{c^2}. \end{equation}

Figure 6.

( $a$ ) Variation of dimensionless recirculation area and lift coefficient with actuation frequency. Panels ( $b$ ) and ( $c$ ) show the spatial distribution of the velocity fluctuation in cases of baseline and optimal open-loop control ( $St=1$ and $D_c=10\,\%$ ).

Figure 6( $a$ ) presents the variations of dimensionless recirculation area $A_{\mathit{b}}^*$ and lift coefficient $C_l$ with Strouhal number. Opposite trends are exhibited by the two parameters: $A_{\mathit{b}}^*$ decreases sharply in the low-frequency range of $0\leqslant St\leqslant 1$ and increases mildly afterwards; $C_l$ , on the contrary, first rises and then drops. The most flow separation elimination is achieved at $St = 1$ , with a minimum value of $A_{\mathit{b}}^* = 0.117$ and a maximum value of $C_l = 0.857$ . Slightly deviating from the optimal Strouhal number (i.e. $0.5 \leqslant St \leqslant 2$ ) would not lead to a drastic decline of the lift coefficient. These observations are consistent with the conclusions in Zong et al. (Reference Zong, Su, Liang and Wu2018), where a completely different actuator (plasma synthetic jet) was used to suppress the leading-edge flow separation of a NACA 0015 aerofoil at $\alpha = 22^\circ$ . Since the flow separation control mechanism of most actuators can be ascribed to enhanced momentum exchange (Greenblatt & Wygnanski Reference Greenblatt and Wygnanski2000; Vernet, Örlü & Alfredsson Reference Tao, Shi, Xiong and Zhou2018; Abe, Nonomura & Fujii Reference Abe, Nonomura and Fujii2023), it becomes natural to also check the turbulent fluctuation in our scenario. Figure 6( $b$ ) compares the spatial distributions of the in-plane velocity fluctuation amplitude (denoted as $\widetilde {U}_{\textit{xy}}$ ) under baseline and optimal open-loop control ( $St = 1$ ). Here, $\widetilde {U}_{\textit{xy}}$ is computed as

(3.4) \begin{equation} \widetilde {U}_{\textit{xy}}=\sqrt {\langle u'u'\rangle +\langle v'v'\rangle } ,\end{equation}

where $\langle u'u'\rangle$ and $\langle v'v'\rangle$ represent the streamwise and wall-normal Reynolds normal stresses, respectively. In both cases, the distribution of $\widetilde {U}_{\textit{xy}}$ takes a conical shape expanding from upstream to downstream, and violent fluctuations are mostly observed in the vicinity of the zero-velocity contour line. Compared with the baseline cases where the peak value of $\widetilde {U}_{\textit{xy}}$ resides at around $(x,\;y)=(0.7c,\, 0.2c)$ , the open-loop control case demonstrates an increased peak velocity fluctuation level ( $0.4\,U_\infty \rightarrow 0.5\, U_\infty$ ) at a more upstream location of $(x,\;y)=(0.5c,\, 0.15c)$ , indicating that the unsteady motions of the shear layer have been excited by the plasma-induced periodical disturbance.

4.1. Implementation of the DDPG algorithm

Unlike the open-loop periodical control in § 3, DRL control requires real-time adjustment of the plasma actuation strength (i.e. action value) according to sensor measurements. For this purpose, the high-speed FPGA-based experimental DRL (FeDRL) control framework proposed in Zong et al. (Reference Zong, van and Kotsonis2025) is employed, and necessary modifications are made to implement a policy-based DRL algorithm (DDPG). The core idea of FeDRL lies in accelerating the inference of neural network control law through an FPGA chip, making it particularly suitable for real-time closed-loop flow control of fast air flows where a control loop frequency of 1–10 kHz is necessitated. The DDPG algorithm is an off-policy reinforcement learning algorithm developed from the actor-critic architecture, capable of solving problems in both continuous and discrete action spaces (Lillicrap et al. Reference Lillicrap, Hunt, Pritzel, Heess, Erez, Tassa, Silver and Wierstra2015). This classical algorithm is preferred, instead of more popular ones such as trust region policy optimisation (Schulman et al. (Reference Roupassov, Nikipelov, Nudnova and Starikovskii2015)) and proximal policy optimisation (PPO, Schulman et al. (Reference Schulman, Levine, Moritz, Jordan and Abbeel2017)), mainly due to the high data efficiency.

Figure 7 sketches the data and work flow in the DRL control. Two parallel loops are configured: a real-time control loop at 1000 Hz and a control law training loop at 3 Hz. Within the control loop, the instantaneous hotwire velocity ( $u_h(t)$ ) is first sampled via an analogue-to-digital converter (ADC), then normalised with the baseline statistics, and finally, concatenated with its historical values to form a state vector ( $\boldsymbol{x}(t)$ ). Mathematically, we have

(4.1) \begin{equation} \left \lbrace \!\!\begin{array}{l} s(t)=\displaystyle \frac {u_h(t)-U_{h,\textit{BL}}}{\sigma _{h,\textit{BL}}}, \\[12pt] \boldsymbol{x}(t)=[s(t), s(t-T_c), s(t-2 T_c), s(t- (n-1)T_c)], \end{array} \right . \end{equation}

where $U_{h,\textit{BL}}$ and $\sigma _{h,\textit{BL}}$ denote respectively the mean and standard deviation of the hotwire velocity measured under baseline conditions; $n$ is the state vector dimension and $T_c$ is the control period ( $1/f_c$ ). Assuming that $u_h(t)$ follows approximately a Gaussian distribution, the dimensionless hotwire velocity $s(t)$ will have a typical range of [−3, 3]. For brevity, $s(t)$ and $\boldsymbol{x}(t)$ are noted as $s_t$ and $\boldsymbol{x}_t$ hereinafter. According to Takens theorem (Takens Reference Suárez, Alcántara-Ávila, Miró, Rabault, Font, Lehmkuhl and Vinuesa2006), the complete dynamic behaviour of high-dimensional systems can be reconstructed from low-dimensional time series. Therefore, augmentation of the state vector with historical measurements is a common practice in experimental ML control (Li et al. Reference Li, Noack, Cordier, Borée and Harambat2017a ; Xia et al. Reference Wang, Mei, Aubry, Chen, Wu and Wu2024), where only a few sensors can be configured in the flow field. While feeding the historical sensor measurements sequentially into a long short-term memory (LSTM) model can be a feasible alternative to state augmentation, it introduces a non-trivial computational overload. Specifically, for the LSTM model, inference of each time step requires four dense matrix multiplications, resulting in a significant control latency prohibitive for real-time flow control.

Figure 7.

Implementation of the DDPG algorithm using the FeDRL control framework.

Figure 7. Long description

A diagram of the implementation of the DDPG algorithm using the FeDRL control framework. The diagram represents a system for controlling flow separation over an aircraft aerofoil. The overall structure includes several key components: Action, Flow environment, State, Power supply, DAC, ADC, FPGA, and CPU. The diagram shows the interaction between these components in a control loop and a train loop. The control loop operates at 100 microseconds, while the train loop operates at 500 milliseconds. The actor network is implemented by RBF, and the critic network is implemented by MLP. The diagram includes labels such as DAC, ADC, FPGA, CPU, and various network components like the actor network and critic network. The flow of data and control signals is depicted with arrows indicating the direction of information flow. The diagram also includes components like the power supply, sensor signal, and experience data, which are part of the overall control and training process.

As a result, the closed-loop control law is parameterised by a RBF neural network, which takes the augmented state vector as the input and outputs a continuous action value (denoted as $a_t$ ). The relationship between $\boldsymbol{x}_t$ and $a_t$ are as follows:

(4.2) \begin{equation} \left \lbrace\!\! \begin{array}{l} a_t=\mathrm{sigmoid} (\boldsymbol{W}_{1\times m}\boldsymbol{\cdot }\left [ \phi _1(\boldsymbol{x}_t),\phi _2(\boldsymbol{x}_t),\ldots ,\phi _m(\boldsymbol{x}_t)\right ]^{\boldsymbol{T}}+\boldsymbol{B} ),\\[3pt] \phi _i=\mathrm{exp}\big(-h_i\left |\left |\boldsymbol{x}_t-\boldsymbol{c_i}\right |\right |^2\big), \hspace {10pt}\mathrm{for} \quad i=[1,2,\ldots ,m]. \end{array} \right . \end{equation}

Here $\phi _i$ is the Gaussian basis function; $\boldsymbol{c_i}$ and $h_i$ are the RBF centre and variance, respectively; $m$ denotes the number of hidden neurones, fixed as 32; $\boldsymbol{W}_{1\times m}$ and $\boldsymbol{B}$ represent the weights and biases connecting the hidden layer to the output. A sigmoid activation function is adopted by the final layer to limit the range of $a_t$ to 0–1. Inspired by Chen et al. (Reference Chen, Zong, Song, Wu, Liang and Xiang2024a ), the output of the actor network $a_t$ is sent to the digital-to-analogue converter (DAC) on the FPGA-based control platform, where the value of $a_t$ is further mapped to the instantaneous plasma actuation strength by pulse width modulation. The DAC converts the digital control signal $a_t$ into an analogue voltage that drives the power supply, and the switch-on time of the enable signal in one control period is given by $t_{en}=a_t\times T_c$ . This control set-up is different from that used in Zong et al. (Reference Zong, van and Kotsonis2025), where a discrete logic control was performed, and the network outputs represent Q values of different actions.

A Gaussian noise (denoted as $\epsilon$ , $\epsilon \sim N(0,\sigma ^2)$ ) is superimposed on the network’s output to explore the unknown state-action space. During the course of learning, the standard deviation of the noise ( $\epsilon$ ) decreases linearly from an initial value of 0.3 towards zero. Such a variation is understandable since at the end of training, the agent becomes more confident about its choice and no more exploration is needed. After each control, the state and action values are organised into an experience tuple $(\boldsymbol{s_t},a_t,\boldsymbol{s_{t+1}})$ and transmitted in batch to the host via a first-in-first-out (FIFO) sequence. Subsequently, the host computes the rewards $r_t$ of that data sequence according to (4.3) and stores the complete control experience to a replay buffer. Finally, mini-batches are randomly sampled from the buffer for network training.

Inspired by Xiang et al. (Reference Xia, Zhang, Kerrigan and Rigas2024) and Chen et al. (Reference Chen, Zong, Wu, Ye, Song and Liang2024b ), the reward is designed as the sliding-average value of the dimensionless state in a time-delayed interrogation window, i.e.

(4.3) \begin{equation} r_t=\displaystyle \frac {1}{R_l}\sum \limits _{i=1}^{R_l}s_{t+i+R_d} \end{equation}

where $R_d$ and $R_l$ denote the reward time delay and reward average length, respectively. For the special case of $R_d=0$ and $R_l=1$ , the instant reward is essentially the dimensionless velocity measured at the next time step, which is reasonable since a positive value of $s_t$ indicates a below-average flow separation area. Addition of the sliding-average operation in (4.3) aims to reduce the measurement uncertainty and, thus, improves the stability of learning. The time-delay operation accounts for the time required for a disturbance to propagate from the upstream actuator to the downstream sensor (Belson et al. Reference Belson, Semeraro, Rowley and Henningson2013).

In DDPG the loss of the critic network is constructed from the Bellman equation (Meyn Reference Meyn2024), i.e.

(4.4) \begin{equation} \left \lbrace\!\! \begin{array}{l} Q_{\textit{target}}=r_t+\gamma\! Q^{\prime }\big(\boldsymbol{x}_{t+1},\mu ^{\prime }\big(\boldsymbol{x}_{t+1};\theta ^{\mu ^{\prime }}\big);\theta ^{Q^{\prime }}\big) ,\\[6pt] L\big(\theta ^Q\big)=\displaystyle \frac {1}{N_b}\sum \limits _{i = 1}^{N_b}\big(Q\big(\boldsymbol{\boldsymbol{x}_{t}},a_t;\theta ^Q\big)-Q_{\textit{target}}\big)^2,\end{array} \right . \end{equation}

where $\gamma$ is the discount factor, set as 0.99; $Q$ and $Q^\prime$ denote the critic and target critic networks, respectively; $\mu$ and $\mu ^\prime$ represent the actor and target actor networks, respectively; $L$ is the loss term; $\theta$ generally refers to the network parameter; $N_b$ is the batch size. As evident, the critic network is optimised in such a way that the predicted Q value is close to the target Q value.

For the actor network, the loss is set as the negative of the mean Q value, such that the objective is to train a policy that maximises the expectation of long-term reward:

(4.5) \begin{equation} L(\theta ^\mu )=-\frac {1}{N_b}\sum \limits _{i = 1}^{N_b}Q(\boldsymbol{x}_t,\mu (\boldsymbol{x}_t;\theta ^\mu );\theta ^Q) \end{equation}

Back-propagating the loss terms in the actor and critic networks, the trainable parameters can be updated. Each time after training, the latest coefficients of the actor network, namely $\theta ^{\mu }$ , are first sent back to the LabVIEW program through a TDMS file and then passed to the FPGA controller via a host-to-target FIFO sequence. To ensure the stability of learning, the target network parameters are updated softly as

(4.6) \begin{equation} \left \lbrace \begin{array}{l} \theta ^{Q^{\prime }}\leftarrow \eta \theta ^Q+(1-\eta )\theta ^{Q^{\prime }},\\ \theta ^{\mu ^{\prime }}\leftarrow \eta \theta ^\mu +(1-\eta )\theta ^{\mu ^{\prime }},\\\end{array} \right . \end{equation}

where $\eta$ is a parameter controlling the target update rate (0.05). The hyperparameters of the DDPG algorithm are summarised in table 2.

Table 2.

Hyper-parameters in DRL control.

4.2. Impact of state vector length

Figure 8 shows the time evolution of several key parameters in a representative DRL control process. The input dimension, sensor position, reward time delay and reward smoothing length are set as $n=4$ , $y_h=0.1\,c$ , $R_d=5$ and $R_l=5$ , respectively. Three repetitive measurements are conducted, and each case lasts for a total test time of 50 s. The first 35 s belong to the training stage, while the last 15 s correspond to the evaluation stage where the optimal control law is executed deterministically without random noise. During the course of training, the instant reward $r_t$ first fluctuates at around 0, then increases rapidly between 10 s $\leqslant t\leqslant$ 20 s and eventually stabilise at a peak reward of approximately 0.2 (see figure 8 $a$ ). Under baseline conditions, $\Delta C_l$ fluctuates around 0 as expected. The lift coefficient exhibits a similar variation trend as the instant reward, whereas the dimensionless flow separation area manifests an opposite trend. The time-averaged values of $\Delta C_l$ and $A_{\mathit{b}}^*$ in the evaluation stage reach 0.148 and 0.117, respectively, comparable to that achieved by the optimal open-loop control (see figure 6 $a$ ). The actuation probability ( $P_a$ ) defined as the mean action value within a short time period is further plotted in figure 8 ( $d$ ). As a result, $a_t$ first goes up from 0.5 to 0.8 during the initial learning stage of $t\leqslant$ 10 s, then drops asymptotically towards 0.13 between 10 s $\leqslant t\leqslant$ 30 s and finally remains invariant. This time evolution of $a_t$ reveals that low-intensity plasma actuation is more effective in flow separation control, consistent with the lift increment map in figure 4(b). Based on the above analysis, the most attractive feature of DRL control probably goes to the high efficiency. Only 50 s are needed to find a closed-loop control law that is on par with the optimal open-loop control, eliminating the tedious parameter tuning process.

Figure 8.

Time evolution of (a) the instant reward, (b) lift coefficient, (c) dimensionless flow separation area, and (d) actuation probability in a representative DRL control process. In each plot the thick line and the shaded area represent the mean value and the standard deviation obtained from three repetitive measurements.

The following three factors contribute to the rapid training of the DRL control strategy. First, a novel high-speed control framework is proposed. Compared with the traditional CPU-based DRL framework, the present study adopts a high-speed FeDRL control framework that can achieve extremely high control frequencies (0.1–10 kHZ), enabling a rapid collection of control experience. Second, data efficient value-based DRL algorithms are deployed. Existing numerical studies favour the usage of on-policy reinforcement learning algorithms such as PPO. Although they are robust in training, the control policy can only be updated by the experience collected by itself, resulting in a low data efficiency. As a comparison, the present study adopts an off-policy algorithm (DQN, DDPG) that allows reuse of all the past experience during training. Third, much less neural network parameters are involved. Unlike the widely used multilayer perceptron (MLP) neural network that has a typical parameter count of 1000–10,000, the present study choose to represent the closed-loop control law with a lightweight RBF neural network with a parameter count of the order of 100. Consequently, much less training steps are requested for neural network convergence.

Zong et al. (Reference Zong, Wu, Li, Su and Liang2024) and Wang et al. (Reference Wang, Hua, Aubry, Chen, Wu and Cui2024) found that by incorporating historical measurements of a single sensor as the state input, a full-state observation of the system can be approximated, leading to an increase in the long-term reward. However, at this point, it is not clear what is the optimal amount of historical data that should be included in the state vector (namely the optimal state input length). Intuitively, the instantaneous velocity trace measured by the hotwire reflects the dynamics of turbulent flow. Within a certain time window, the two measured velocity points are correlated, revealing the motion of large-scale structures. In contrast, when two velocity points are obtained too far away from each other in the time frame, no correlation should be expected, and thus, incorporation of more distant sensor information can be useless. Redundant information is regarded as turbulent noise, which may lead to an inaccurate estimation of the system state and a biased decision in closed-loop control (Peng et al. Reference Passaggia, Quansah, Mazellier, Maceda and Kourta2025).

Figure 9.

(a) Time trace of the instantaneous hotwire velocity under baseline conditions. (b) Auto-correlation of the hotwire velocity signal. (c) Increment of the lift coefficient obtained by DRL control at an increasing state input length.

Figure 9. Long description

The image contains three separate graphs. Panel A: A line graph showing the time trace of the instantaneous hotwire velocity under baseline conditions. The x-axis is labeled t/T0 and ranges from 0 to 100. The y-axis is labeled uh/U∞ and ranges from 0 to 1. The graph displays fluctuating velocity values over time. Panel B: A line graph showing the auto-correlation of the hotwire velocity signal. The x-axis is labeled τ/T0 and ranges from 0 to 3. The y-axis is labeled RU and ranges from −0.4 to 1. The graph shows a decreasing trend in auto-correlation over time. Panel C: A line graph showing the increment of the lift coefficient obtained by DRL control at an increasing state input length. The x-axis is labeled n and ranges from 0 to 25. The y-axis is labeled ΔCL and ranges from 0 to 0.15. The graph displays fluctuating values of the lift coefficient increment with error bars indicating variability.

To validate the above conjecture, the temporal evolution of the instantaneous hotwire velocity signal measured under baseline conditions is plotted in figure 9(a). The characteristic time scale of the flow ( $T_0=c/U_\infty$ ) and the control period ( $T_c$ ) are used as the reference scales for time normalisation. The instantaneous hotwire velocity $u_h$ oscillates violently around a mean value of 0.35 $U_\infty$ , and high-rise spikes can be picked occasionally. The oscillation period of $u_h$ seems to differ significantly from cycle to cycle, indicating a broad distribution of the frequency spectrum, consistent with the observation in Romeos et al. (Reference Romdhana, Merlo, Ceccato and Tonella2009). Figure 9(b) shows the auto-correlation coefficient of the hotwire velocity signal ( $R_u$ ), computed using

(4.7) \begin{equation} R_u(\tau )=\frac {1}{N_s\times \sigma _{h,\textit{BL}}}\sum \limits _{t = \tau }^{T_{\textit{max}}}(u_h(t-\tau )-U_{h,\textit{BL}})(u_h(t)-U_{h,\textit{BL}}), \end{equation}

where $\tau$ and $T_{\textit{max}}$ represent the time delay and the total sample time, respectively; $N_s$ is the data sequence length.

In case of zero time delay, $R_u$ equals to 1, by definition. With increasing time delay from 0 to $3T_c$ , the auto-correlation coefficient first drops, then goes up and finally decreases again. After every oscillation, the peak amplitude of $R_u$ decays noticeably. This behaviour is very similar to the sinc function (i.e. $f(t)=sin(t)/t$ ) frequently used in signal processing, indicating that the raw hotwire velocity does exhibit an approximate periodical variation. Based on the positions of local minimum and maximum ( $1.2T_c$ and $2.4T_c$ ), the average oscillation period of the hotwire velocity is estimated to be $2.4T_c$ , corresponding to a characteristic frequency of $0.4f_c$ . According to Chang et al. (Reference Chang, Zhang, He and Zhou2022), the vortex shedding frequency during flow separation of a NACA 0012 aerofoil at different angles of attack is 0.2–0.6 $f_c$ . Such a coincidence implies that the measured periodicity is mostly caused by vortex shedding. Additionally, by integrating the auto-correlation coefficient from $\tau =0$ to the first zero-crossing point of $\tau =8T_c$ , the Taylor microscale of the separated aerofoil flow is determined to be $4T_c$ (Pope Reference Pino, Schena, Rabault and Mendez2000).

Figure 9(c) shows the lift coefficient increment obtained by DRL control at different state vector lengths. In the range of $1\leqslant n\leqslant 4$ , $\Delta C_l$ increases monotonically with the number of state inputs, proving that incorporating appropriate historical information to the flow does improve the DRL performance, consistent with the findings of Wang et al. (Reference Wang, Hua, Aubry, Chen, Wu and Cui2024). However, when further increasing the state vector length from 4 to 12, a drastic degradation of the DRL control performance is observed, and the value of $\Delta C_l$ in the case of $n=12$ (0.03) is approximately a quarter of that achieved at $n=4$ (0.12). Connecting this curve with the auto-correlation coefficient shown in figure 9(b), we observe that the optimal state vector length is essentially equal to the ratio of the Taylor microscale to the control period (i.e. 4). Moreover, the worst value of $n$ approximately corresponds to the lowest point in the auto-correlation curve ( $\tau =12T_c$ ). These connections indicate that inclusion of historical flow information into the state vector is only helpful when the add-on is positively correlated with the present state. Due to the nature of the RBF layer, there is no way to selectively disable the impacts of part of the input vector. Consequently, distant historical measurements that have no relationship with the present state can be treated as random noise, which will jeopardise the performance of DRL by polluting the state observation. In future studies, which involve constructing the actor network, a fully connected layer might be added in between the input vector and the RBF layer. This extra layer is expected to assign the neural network an ability to discard irrelevant input elements during training, thereby eliminating their adverse influence on the network’s output.

Figure 10.

Trajectories of DRL control at different input lengths. From left to right, the values of $n$ in each column are 1, 4 and 12. From top to bottom, the plotted variables in each row are the instant reward, lift coefficient increment and actuation probability. The same plotting style as that in figure 8 is used.

Figure 10. Long description

The image contains nine line graphs arranged in a 3x3 grid. Each column represents a different input length value (1, 4, and 12 from left to right). Each row represents different plotted variables: instant reward, lift coefficient increment, and actuation probability from top to bottom. Panel A: The first row shows the instant reward (r_t) over time (t in seconds). The y-axis ranges from −0.5 to 0.5. Panel B: The second row shows the lift coefficient increment (ΔC_l) over time (t in seconds). The y-axis ranges from 0 to 0.20. Panel C: The third row shows the actuation probability (P_a) over time (t in seconds). The y-axis ranges from 0 to 1.0. Each graph includes a line with a shaded area representing variability or confidence intervals.

The trajectories of DRL control at three representative state vector lengths of $n = 1, 4, 12$ are further compared in figure 10. When the state vector contains only one element, the instantaneous reward, lift coefficient and actuation probability exhibit significant fluctuations, rendering the entire training process unstable. This phenomenon was also reported in the Q-learning control of Zong et al. (Reference Zong, Wu, Li, Su and Liang2024), where they failed to get a converged control law when using only the present velocity as the state vector. As a comparison, in the optimal case of $n=4$ , due to the inclusion of correlated historical flow information, both $r_t$ and $\Delta C_l$ increase steadily. The lift coefficient increment in the evaluation stage reaches 0.12 at a low-actuation probability of $P_a=0.21$ . When the value of $n$ is boosted to 12, uncorrelated redundant information starts to be blended in the sate vector. Consequently, the DRL agent fails to pick out the states that require plasma actuation, and the mean actuation probability decreases linearly towards zero during the training process. The lift coefficient increment earned by DRL is only 0.03.

4.3. Impact of sensor placement

In the DRL-based flow control, sensor configuration is crucial, as the premise of optimal control is the awareness of dynamics in the flow (Xu & Zhang Reference Xiang, Zong, Wu, Li and Liang2023). For numerical studies, near-complete observation of the flow system can be achieved by deploying hundreds of sensors in the flow fields (Rabault et al. Reference Pope2019; Wang et al. Reference Wang, Song, Cai, Liu and Wen2022b ; Chen et al. Reference Chen, Zong, Song, Wu, Liang and Su2023). As a comparison, in experimental studies, only a few feedback sensors can be configured in the flow (Amico et al. Reference Amico, Cafiero and Iuso2022; Zong et al. Reference Zong, van and Kotsonis2025), and the position of the feedback sensor really impacts the performance of DRL.

Figure 11.

( $a$ , $b$ ) Contours of the mean velocity and velocity fluctuation measured by PIV under baseline conditions. Purple dots indicate the sensor positions. ( $c$ , $d$ ) Vertical profiles of the mean velocity and velocity fluctuation extracted at $x=1.1c$ . The purple line shows the lift coefficient increment achieved by DRL at different sensor positions.

Huerre & Monkewitz (Reference Huerre and Monkewitz1990) proposed that the flow dynamics system fluctuates at a specific frequency, and the major disturbances are always located in a specific region of the flow field. Figures 11(a) and 11(b) plot the mean velocity and fluctuation velocity contours measured by PIV under baseline conditions. As indicated by the pink dots, six different sensor positions are realised with the traversing system, covering the entire shear layer. The $x$ coordinates of these positions are fixed at $x=1.1\,c$ , and the $y$ coordinates are distributed uniformly within $-0.1\,c\leqslant y\leqslant 0.23\,c$ with a spacing of $0.067\,c$ (10 mm). At each position, an independent DRL control is conducted. For the innermost position of $y=-0.1\,c$ , the hotwire sensor is fully submerged in the recirculation zone. As the hotwire moves upwards, it first probes the separated shear layer and then dives into the free stream. For reference, the mean velocity and fluctuation velocity profiles at $x=1.1c$ are extracted from the contours and further shown in figures 11(c) and 11(d). In the range of $-0.1\,c\leqslant y\leqslant 0.23\,c$ , the mean velocity first drops and then goes up rapidly, with the minimum value achieved at the centre of the shear layer (approximately $y=-0.05\,c$ ). Correspondingly, the velocity fluctuation exhibits a mild increase followed by a sharp decrease, which is reasonable since the majority of the velocity fluctuations are caused by the shedding of shear layer vortices.

The lift coefficient increments obtained by DRL at different sensor positions are plotted side by side with the velocity fluctuation in figure 11( $d$ ). Obviously, there is a positive correlation between $\widetilde {U}_{\textit{xy}}$ and the value of $\Delta C_l$ . The optimal sensor position in DRL coincides with the point of highest velocity fluctuation ( $y=0.03\,c$ ). From the perspective of informatics, fluctuation means dynamics of the controlled system, which is necessary for the agent to understand the behaviour of the flow. In the study by Parezanović et al. (Reference Özkan2016), the hotwire array was placed in the mixing layer to extract the major fluctuation information, and the lock-on state was excited through a closed-loop control based on GP. Considering the extreme case of placing the senor in the free stream, the dynamics of flow separation is completely missing and no effective control strategy shall be expected from the DRL agent.

Figure 12.

Trajectories of DRL control at different sensor placements. From left to right, the positions of the hotwire in each column are −0.1 $c$ , $0.1c$ and $0.23c$ . The same plotting style as that in figure 8 is used.

Figure 12. Long description

The image contains nine line graphs arranged in a 3x3 grid, each representing trajectories of DRL control at different sensor placements. Panel A, Panel B, and Panel C show data for different positions of the hotwire. Each panel contains three subplots labeled (a), (b), and (c) for the top row, (d), (e), and (f) for the middle row, and (g), (h), and (i) for the bottom row. The x-axis represents time (t) in seconds, and the y-axis represents different variables: r_t for the top row, ΔC_l for the middle row, and P_a for the bottom row. The blue, red, and magenta lines indicate different data series. The shaded regions around the lines represent the uncertainty or variability in the data. Each subplot shows trends over time, with varying patterns and fluctuations. The graphs illustrate how the trajectories change with different sensor placements.

Figure 12 compares the time histories of DRL control at three representative sensor locations. When the hotwire is placed in the recirculation zone ( $y = -0.1\,c$ ), unstable learning with opposite variations of $r_t$ and $\Delta C_l$ are manifested, unlike what we have observed in figure 10. Specifically, during the initial training stage, trail operation of the plasma actuator increases the lift coefficient by a considerable amount, achieving flow separation suppression. Such a shrinkage of the flow separation area, however, reduces the velocity measured by the hotwire, leading to a negative reward because the separated shear layer gets closer to the aerofoil suction surface (see figure 5 $d$ ). During DRL, the control strategy is updated based on the reward-maximisation principle. Therefore, the initial ‘successful’ trail operations are mistakenly labelled as harmful actions by the agent, and the actuation probability is continuously cut down to zero until the end of training. In the case of sensor placement at the free stream ( $y = 0.23\,c$ ), similar variation trends of the instant reward and lift coefficients are exhibited, which once again illustrates that the premise of successful DRL control is the sensation of flow dynamics. As a comparison to the above two cases, an optimal sensor placement of $y = 0.1\,c$ gives a steady and concurrent rise of $r_t$ and $\Delta C_l$ .

Based on the above observation, it can be concluded that for simple two-dimensional flows like aerofoil flow, an appropriate layout of one sensor is already sufficient to perform successful closed-loop control. A similar practice can be found in Maceda et al. (Reference Maceda, Varon, Lusseyran and Noack2023) where only one hotwire is placed in the shear layer of the open cavity, acting as the feedback sensor for the LGP control. However, for complex three-dimensional flows (such as hypersonic inlet, Ahmed body or multi-element wing), a sophisticated feedback sensor configuration is required to fully capture the instantaneous state of the flow (Amico et al. Reference Amico, Cafiero and Iuso2022; Tao et al. Reference Takens2025).

4.4. Impact of reward shaping

In this study, since the disturbances induced by an upstream plasma actuator require some time to propagate to the downstream hotwire sensor, a time delay inevitably exists between the action and its corresponding reward. For non-real-time flow control, this time delay can be largely overlooked since no instant sensor feedback is required for the action. Whereas, for real-time closed-loop control, ignoring the delay effect will be detrimental to the control performance and affect the stability of the control system (Zhu & Krstic Reference Zhou, Xia, Luo, Wang and Deng2020).

Figure 13.

(a) Time evolution of the instantaneous velocity and actuation signal in a threshold control test. (b) Impulse responses of the instantaneous hotwire velocity to plasma actuation collected during 100 repetitive measurements (thin grey lines). The thick black line represents the ensemble-averaged result.

Brzozowski et al. (Reference Brzozowski, Woo, Culp and Glezer2010) found that instantaneous pulse actuation can lead to significant transient change in flow dynamics, and the changes can be quantitatively evaluated through phase-locked measurements. In this study, to quantify the time delay, the impulse response of the hotwire velocity to a single shot of plasma actuation is acquired by performing a threshold control test. Figure 13(a) shows representative waveforms of the hotwire sensor signal and the actuation command during such a test. Specifically, when the dimensionless hotwire velocity $s_t$ positively crosses a threshold value (denoted as $S_{th}$ , 1.5 in this example), the actuator immediately shoots a pulse (width: 1 ms). Then, the impulse response of the flow is observed quietly for 100 ms (10 $T_0$ ) before the next trigger pulse can be activated. Figure 13(b) shows the instantaneous hotwire sensor signals collected during 100 threshold control tests. At the trigger moment of $t/T_0=0$ , the hotwire signals exhibit low variance since the threshold condition has to be satisfied. Moving towards the post-actuation time zone, significant velocity fluctuations are observed, ranging from -2.5 to 2.5. To extract the mean variation trend, an ensemble-averaging operation is applied to these curves, resulting in a time-averaged impulse response curve shown in thick black. As expected, the mean hotwire velocity gradually recovers from the threshold level of 1.5–0.

Figure 14(a) displays the impulse response curves obtained at decreasing threshold control levels of $S_{th}=0.5,0,-0.5\,\mathrm{and} -1$ . In each plot, the natural evolution trend of the hotwire velocity without plasma actuation is plotted side by side for comparison. For all the cases, a typical time delay between 0.2 $T_c$ and 0.5 $T_c$ is required for the effect of plasma actuation to be discerned, justifying our previous discussion of the reward delay. Moreover, the impact of a single-shot plasma actuation can last for more than 3 $T_c$ . By subtracting the baseline curve from the actuated one, the differential velocity is found to be first positive ( $0.2\lt t/T_c\leqslant 1.2$ ), then negative ( $1.2\lt t/T_c\leqslant 1.8$ ) and finally positive again ( $1.8\lt t/T_c\leqslant 3.0$ ), indicating an undulating motion of the separated shear layer. For stable training in DRL, these two important features of the reward, namely ‘significant time delay’ and ‘long cashing period’, must be taken into consideration (Chen et al. Reference Chen, Zong, Wu, Ye, Song and Liang2024b ).

Figure 14.

(a) Impulse response curves of the hotwire velocity at decreasing threshold control levels of 0.5, 0, −0.5 and −1. (b) The lift coefficient increment map spanned by the reward time delay and reward smoothing length. Purple stars indicate the locations of test points.

Figure 14. Long description

Panel A: Four line graphs show impulse response curves of the hotwire velocity at decreasing threshold control levels of 0.5, 0, −0.5, and −1. The x-axis represents normalized time (t/T0) ranging from 0 to 3, and the y-axis represents the velocity (st) ranging from −1 to 1. The red squares represent threshold control, and the blue circles represent baseline. Each graph shows different trends in velocity response over time for each control level. Panel B: A contour plot displays the lift coefficient increment map spanned by the reward time delay (Rd*) on the y-axis ranging from 0 to 2 and the reward smoothing length (Rl*) on the x-axis ranging from 0 to 2. The color gradient from blue to red indicates the lift coefficient increment (ΔCl) values, with blue representing lower values and red representing higher values. Purple stars mark the locations of test points across the map. The contour lines highlight regions of similar lift coefficient increments.

Figure 15.

Trajectories of ( $a$ ) the instant reward, ( $b$ ) lift coefficient increment, and ( $c$ ) actuation probability in DRL control under a typical reward setting of $R_d^\ast =1.5$ and $R_l^\ast =0.5$ . The same plotting style as that in figure 8 is used.

The setting of instantaneous rewards is an important factor restricting the application of DRL control. For complex flow scenarios, improper reward setting may lead to ineffective control strategies (Vignon, Rabault & Vinuesa Reference Vernet, Örlü and Alfredsson2023). To examine the impact of reward settings, the two parameters in (4.3) ( $R_d$ and $R_l$ ) are taken as decision variables. For each variable, five discrete values are selected, resulting in a total of 25 parameter combinations. For each combination, three repeated DRL controls are conducted, and the average $\Delta C_l$ during the evaluation stage is plotted in figure 14(b). Since the reward time delay is in units of control period, $T_c$ and $T_d$ can be used to normalise $R_d$ , resulting in a dimensionless reward time delay of $R_d^\ast =R_dT_c/T_0$ . By analogy, a dimensionless reward smoothing length is derived from $R_l^\ast =R_lT_c/T_0$ . Based on the lift increment map, the reward time delay plays a more influential role than the reward smoothing length. Particularly, at a combination of $R_d^\ast = 0.5$ and $R_l^\ast = 0.5$ , DRL control achieves the most lift coefficient increment (0.148). The instant reward at this time essentially corresponds to the average value of $s_t$ between $0.5\lt t/T_c\lt 1$ (see figure 14 a), leading to a positive reward of plasma actuation. Consequently, a stable learning process is exhibited, as illustrated in figure 8. For the case of $R_l^\ast = 0.5$ and $R_d^\ast = 1.5$ , DRL control only manages to increase the lift coefficient by 0.08. This is not surprising, since the average value of $s_t$ between $1.5\lt t/T_c\lt 2$ is taken as the instant reward, which according to figure 14( $a$ ) assigns approximately zero reward to plasma actuation. The corresponding history of DRL control is plotted in figure 15. Consequently, the agent can hardly achieve any positive rewards, and the actuation strength undergoes significant fluctuations, indicating an unsuccessful learning that failed to get any converged control laws.

Despite the considerable efforts made to optimise the state input, sensor position and reward function, the performance of the optimal control law is still not able to surpass the best open-loop control. From the authors’ perspective, whether a closed-loop control strategy that surpasses the best open-loop control strategy can be discovered depends more on the flow that is controlled, instead of the ML control algorithm that is adopted. For cases where actuation should be applied conditionally on the phase/event of the flow, closed-loop control is expected to win out. A typical example is open cavity flow (Yan et al. Reference Yagiz, Kandil and Pehlivanoglu2006), where manually designed closed-loop control laws as simple as proportional-integral-derivative can already outperform largely the best open-loop control, not to mention advanced control algorithms such as LGP and DRL. In contrast, for cases where periodical actuation is sufficient to excite the flow (e.g. leading-edge aerofoil flow separation and cylinder drag reduction), the optimal closed-loop control laws discovered by ML methods can only deliver comparable performance as the best open-loop control (Fan et al. Reference Fan, Yang, Wang, Triantafyllou and Karniadakis2020; Shimomura et al. Reference Schulman, Wolski, Dhariwal, Radford and Klimov2020).

5.1. Low input dimension of $n\leqslant 2$

Figure 16.

(a) The state distribution probability and the final neural network control law obtained by DRL at $n=1$ . (b) Performance comparison between the neural network control law and the polynomial control law. The dash-dotted red line indicates the control start moment.

Figure 16. Long description

The image contains two line graphs labeled as Panel A and Panel B. Panel A shows a comparison between neural network output, polynomial output, and state probability. The horizontal axis is labeled with s_t and ranges from −3 to 3. The left vertical axis is labeled a_t and ranges from 0 to 0.5, while the right vertical axis is labeled P and also ranges from 0 to 0.5. The neural network output is represented by blue circles connected by a blue line, the polynomial output by a solid blue line, and the state probability by a dashed magenta line. The neural network and polynomial outputs show a sinusoidal pattern, while the state probability forms a bell-shaped curve. Panel B compares the performance of the neural network control law and the polynomial control law over time. The horizontal axis is labeled t (s) and ranges from 0 to 30 seconds. The vertical axis is labeled ΔC_l and ranges from −0.10 to 0.20. The neural network control law is represented by a blue line with a shaded blue area indicating variability, and the polynomial control law by a red line with a shaded red area. Both lines show an increase in performance after the control start moment, indicated by a black dash-dotted vertical line at around 15 seconds.

Interpretation of the control strategy essentially requires us to understand how the neural network maps a state vector to the control action ( $\mu _\theta :\boldsymbol{x}_t\rightarrow a_t$ ). For a low-dimensional state input of $n\leqslant 2$ , the relationship between $\boldsymbol{x}_t$ and $a_t$ can be directly visualised by a curve or contour. Figure 16(a) plots the final neural network control law obtained by DRL at $n=1$ . In this case, the state vector $\boldsymbol{x}_t$ incorporates only one element of $s_t$ , and the probability distribution of $s_t$ resembles very much a standard normal distribution, with the majority of the states residing in the range of $-2\leqslant s_t\leqslant 2$ . The action value $a_t$ seems to be a symmetric function of the state variable, and the actuation intensity is approximately 0 between $-1\leqslant s_t\leqslant 1$ . As the absolute value of $s_t$ goes beyond 1, the action value starts to rise rapidly, reaching a maximum value of 0.43. Since $s_t$ is the normalised hotwire velocity reflecting how much the instantaneous velocity deviates from the baseline mean, the one-dimensional (1-D) control law dictates that strong plasma actuation should only be imposed when a severe expansion/contraction of the flow separation zone occurs. This strategy parallels the closed-loop approach of Gautier et al. (Reference Gautier, Aider, Duriez, Noack, Segond and Abel2015), who used the separation zone area measured by PIV as a feedback signal to control flow separation in the water flow.

The control logic at $n=1$ is so simple that it is possible to distil the neural network control law into a polynomial function. Such a simplification can not only reduce the computation cost, but also improve the maximum control loop frequency. As shown in figure 16(a), a 16th-order polynomial is sufficient to accurately reproduce the DRL control law, and the parameter usage is only 20 % of the original RBF network. To validate whether such a simplification can retain the control performance, a step control test is performed, and the time evolution of the lift coefficient increment is compared in figure 16(b). During the first 15 s, control is not applied, and $\Delta C_l$ exhibits mild oscillations around 0. Upon plasma actuation, both control laws produce nearly identical responses of $\Delta C_l$ , characterised by an sharp transient jump followed by a plateau.

Figure 17.

(a) Probability density distribution of the state vector in a two-dimensional space. (b) Visualisation of the two-dimensional neural network control law. (c) Polynominal fitting of the neural network control law. (d) Performance comparison between the neural network control law and the fitted polynomial control law. The red dotted line indicates the control start moment.

In the case of $n=2$ , the probability density distribution of the state vector ( $\boldsymbol{x}_t=[s_t,\,s_{t-1}]$ ) is shown as a contour in figure 17(a). As a result, an approximately elliptical distribution of the state vector is exhibited, with the diagonal line of $s_t=s_{t-1}$ severing the major axis. This particular distribution indicates that the dimensionless velocities acquired at the present and previous moments (i.e. $s_t$ and $s_{t-1}$ ) are largely correlated, which is reasonable since the control period is one order lower than the characteristic time scale of the flow and the measured velocity can only vary marginally in one control period. The final neural network control law obtained by DRL at $n=2$ is visualised in figure 17(b). Similar to the 1-D case, plasma actuation is only imposed when $\boldsymbol{x}_t$ is far away from the origin, and the strength of the plasma actuation goes up rapidly when the state vector norm $||\boldsymbol{x}_t||$ exceeds a threshold value of 2.3. A subtle distinction from the 1-D control law manifests when $n=2$ , the agent prefers to generate stronger actuation when the separation region contracts strongly. Quantitatively, the control amplitude at $s_t=s_{t-1}=2$ exceeds that observed at $s_t=s_{t-1}=-2$ . Such a distribution can be well approximated by the following second-order formula:

(5.1) \begin{equation} a_t=\mathrm{{max}}\big\{0.023*\big[(s_t+0.24)^2+(s_{t-1}+0.2)^2\big]-0.1, 0\big\} .\end{equation}

Here the four coefficients are fitted by the SciPy package in Python. As illustrated in figures 17(c) and 17(d), the outputs of the formula agree well with the neural network predictions, and their control performances are equivalent in terms of the final lift coefficient increment.

5.2. High input dimension of $n=4$

Figure 18.

(a) State vectors represented by each cluster centroid. (b) Proximity map of the state vector. (c) Comparison of the cluster distribution probability in baseline and DRL control.

Figure 18. Long description

Panel A: Ten line graphs show state vectors represented by each cluster centroid over time. Each graph has time (t) in milliseconds on the horizontal axis and state vector (s_t) on the vertical axis. The graphs are labeled with cluster numbers from 1 to 10, each with a distinct color. Panel B: A scatter plot displays the proximity map of the state vector with two axes labeled gamma 1 and gamma 2. Each cluster is represented by a distinct color and labeled with cluster numbers from 1 to 10. Panel C: A bar graph compares the cluster distribution probability in baseline and DRL control. The horizontal axis is labeled Cluster with numbers from 1 to 10, and the vertical axis is labeled P_s. The graph includes two sets of bars for each cluster: green bars for baseline probability (P_s,baseline) and blue bars for actuation probability (P_s,actuation).

Figure 19.

(a,b) Cluster transition probability matrices in baseline and DRL control. (c) Comparison of the cluster-specific reward distributions.

Although insights have been gained by visualising the neural network control law obtained at low input dimensions, interpretation of the final control law obtained by DRL at a high input dimension of $n=4$ remains challenging (see § 4.2). In Li et al. (Reference Li, Zong, Wu, Su and Liang2021) and Maceda et al. (Reference Maceda, Varon, Lusseyran and Noack2023), cluster-based methods were adopted to analyse the close-loop control results of laminar mixing layer and cavity flow, providing an intuitive explanation of the complex system dynamics. Inspired by their studies, a sample set $S$ is first constructed, including the state vectors obtained under both baseline and DRL control ( $n=4$ ). The set $S$ is then partitioned into 10 distinct clusters using the $k$ -means algorithm (Lloyd Reference Lloyd1982). Figure 18(a) shows the state vectors represented by each cluster centroid. According to (4.1), clusters 1–3 can be understood as the contracted state of the separation region, clusters 5–7 as the expanded state of the separation region and the other states are transitional states between the two. Subsequently, both the state vectors and the cluster centroids are projected onto a two-dimensional feature space spanned by the principal components ( $\gamma _1$ and $\gamma _2$ ), via multi-dimensional scaling. This dimensionality reduction preserves the essential topological relationships, meanwhile enabling effective visualisation of the high-dimensional control landscape (Li et al. Reference Li, Tan, Gao and Noack2022). As shown in figure 18(b), the clusters are distributed in an oval shape, pivoting the centre cluster 9. This spatial distribution pattern suggests potential periodicity in the system dynamics, as will be elaborated in the following content (Castellanos et al. Reference Castellanos, Cornejo Maceda, De La Fuente, Noack, Ianiro and Discetti2022). Figure 18(c) plots the cluster distribution probabilities in both baseline and DRL control. Compared with the former, the occurrence probabilities of clusters 1–3 under plasma actuation are elevated noticeably, while those of clusters 4–10 are reduced moderately.

The variation of cluster distribution can be traced back to the change of cluster transition. Denoting $p_{i,j}$ as the probability of a state vector transitioning from cluster $i$ to cluster $j$ , a two-dimensional matrix of $p_{i,j}$ can be constructed, which directly reflects the cluster transition dynamics. Figure 19(a) shows the cluster transition probability matrix in the baseline case. The values of $p_{i,j}$ are directly labelled in the grid, except for the inactive transitions with $p_{i,j}$ less than 0.05. As evident, columns 3, 4, 8, 9, 10 have at least four active elements, indicating that the corresponding clusters are the frequently visited destinations. No wonder that in figure 18( $c$ ) a higher distribution probability is exhibited in clusters 3, 4, 8, 9, 10. Under DRL control (figure 19 $b$ ), although the number of active elements in columns 1–3 does not show obvious increase, the self-transition probability in clusters 1–3 are elevated unexceptionally, which explains the variation of distribution probability after plasma actuation. Although the cluster transfer dynamics have been modified considerably by DRL control, it remains unclear to us how these modifications can lead to a higher reward. To resolve this puzzle, the distributions of instant reward at different clusters are plotted in figure 19( $c$ ). Obviously, clusters 1–3 are the favourable clusters receiving high positive rewards, whereas clusters 5–7 are the unfavourable clusters exhibiting low negative rewards. Regardless of plasma actuation, clusters with positive mean rewards tend to have a narrower distribution of $r_t$ than those with negative mean rewards. The reward-earning mechanisms of DRL control are two fold. On one hand, by modifying the instant reward distribution, the mean reward pertaining to the majority of the clusters (clusters 1, 3, 4, 8–10) are elevated. On the other hand, by intervening the cluster transfer dynamics, the occurrence likelihoods of the high-reward clusters (clusters 1–3) are increased.

Figure 20.

(a) Sketch of the cluster transfer graph. Beside each cluster the left bar indicates the relative magnitude of the reward with respect to the global mean reward (negative: blue, positive: red). The right bar represents the relative actuation probability with respect to the global mean actuation probability (negative: green, positive: pink). The thick black arrows and the thin grey arrows indicate the primary transfer routes with $p_{ij}\geqslant 0.3$ and the secondary transfer routes with $0.1\leqslant p_{ij}\lt 0.3$ , respectively. (b,c) Plasma-induced spanwise velocity difference contours pertaining to clusters 1 and 7. The solid black line represent the zero-velocity lines.

Figure 20. Long description

Panel A: A cluster transfer graph showing various clusters connected by arrows. Each cluster is represented by a colored circle and has two bars beside it. The left bar indicates the relative magnitude of the reward with respect to the global mean reward, with blue representing negative values and red representing positive values. The right bar represents the relative actuation probability with respect to the global mean actuation probability, with green representing negative values and pink representing positive values. Thick black arrows indicate primary transfer routes, while thin grey arrows indicate secondary transfer routes. Panel B: A contour plot showing the plasma-induced spanwise velocity difference for cluster 1. The solid black line represents the zero-velocity line. Panel C: A contour plot showing the plasma-induced spanwise velocity difference for cluster 7. The solid black line represents the zero-velocity line.

Integrating the above analyses, the flow dynamics under DRL control can be summarised into a directional graph model, as shown in figure 20(a). For each cluster, the relative magnitudes of the reward and actuation probability (denoted as $s_i$ and $a_i$ ) are plotted as two bars. The entire graph can be roughly separated into two distinct zones: a high-actuation high-reward zone on the left and a low-actuation low-reward zone on the right. In the left zone (pink shade), the clusters are tightly spaced, and clusters 1–3 form a closed-loop cycle with self-transition as the primary route, suggesting that DRL control strives to trap the flow states within the high-reward loop by sustained strong plasma actuation. For the right zone (green shade), below-average rewards are exhibited in clusters 5–7, which correspond to the scenario of large flow separation. Under such circumstances, the DRL agent chose to apply minimal control to the flow. This strategy is different from the low-dimensional control laws explained in § 5.1 where strong plasma actuation is also imposed at low values of $s_t$ (low reward). Presumably, the DRL agent wants to leverage the natural transfer route of the flow for recovery ( $5\rightarrow 9\rightarrow 10$ or $6\rightarrow 7\rightarrow 8\rightarrow 10$ ).

The above-mentioned closed-loop control logic might be generalised to aerofoil flow separation at different wind speeds, as long as a similar topology is presented. Such generalisation evidence can be found in Parezanović et al. (Reference Özkan2016), where ML control utilising LGP was conducted to promote the mixing of a free shear layer. Their results demonstrated that when the free-stream velocity is slightly increased, the effectiveness of the pre-trained control law can still maintain. However, when the flow condition undergoes significant changes, such as varying angles of attack, the pre-trained control strategy is likely to fail. In this scenario, not only fine tuning but also relocation of the sensor is needed. Limited by the scope of the present paper, such validation experiments are not conducted.

Since PIV measurements are also performed in DRL control, a cluster-averaged velocity field can be readily derived from the raw snapshots and the phase information. Subtracting the baseline velocity field from the cluster-averaged velocity field, we can get a differential velocity field reflecting the impact of plasma actuation. Denoting $\Delta V$ as the vertical component of the differential velocity between the time-averaged flow and the cluster-averaged flow, we have the relation

(5.2) \begin{equation} \Delta V=V_i-V ,\end{equation}

where $V_i$ represents the vertical velocity component in the i th cluster-averaged flow field. Figures 20(b) and 20(c) shows the contours of $\Delta V$ pertaining to clusters 1 and 7, respectively. As evident, alternating patches of negative and positive values of $\Delta V$ are observed in both cases, corresponding to the periodical vortex shedding. The high reward in cluster 1 eventually results in a much smaller flow separation area than cluster 7.

6.1. Control effectiveness

A side-by-side comparison between the valued-based DQN algorithm and the policy-based DDPG algorithm is first conducted, since they belong to the same category of DRL. Generally speaking, DQN is easy to implement and can only deal with discrete control tasks, whereas DDPG is more robust and capable of performing both discrete and continuous controls. Successful application of DQN in flow separation control can be found in Shimomura et al. (Reference Schulman, Wolski, Dhariwal, Radford and Klimov2020) and Xiang et al. (Reference Xia, Zhang, Kerrigan and Rigas2024). To compare theses two algorithms fairly, the same state vector length ( $n=4$ ), sensor placement ( $y_h=0.1\,c$ ), reward setting ( $R_d=5$ , $R_l=5$ ) and inflow conditions are used. Detailed implementation of the DQN can be found in our previous works (Chen et al. Reference Chen, Zong, Wu, Ye, Song and Liang2024b ; Zong et al. Reference Zong, van and Kotsonis2025).

Figure 21 shows the time evolution of $\Delta C_l$ and $P_a$ in DQN and DDPG control. Under DQN control, a converged control law is reached after only 10 s of training, much shorter than that consumed by DDPG (20 s). Such a fast optimisation speed is largely due to the fact that the control task in DQN is relatively easy, since only two choices need to be selected. In contrast, a seamless tuning of the action value is required in DDPG. The final control benefit earned by DQN ( $\Delta C_l=0.148$ ) is comparable to that in DDPG, yet the mean actuation probability (0.4) is almost three times higher than the later (0.13). Since the mean actuation probability directly determines the actuator power consumption, the control strategy in DDPG is obviously more efficient, capable of applying the right amount of plasma force without sacrificing the control performance. Physically, DQN can only impose discrete on/off actuation commands, which can be suboptimal, since not all actuated flow scenarios require a full-throttle plasma force. The DDPG eliminates this weakness by applying the right amount of force, which considerably cuts down the energy consumption. To summarise, in the scenario of aerofoil flow separation control, twice training time is needed for the DDPG to find a more budget control law with a comparable lift increment as the DQN.

Figure 21.

Time evolution of (a) the lift coefficient increment and (b) actuation probability in DQN and DDPG control.

Figure 22.

The evolution of (a) the cost function, (b) the lift coefficient increment, (c) the actuation probability, and (d) the individual control laws during LGP. The symbol $i$ on the $x$ axis represents the number of evaluations. Each point represents an individual control law. In panel (d) the points are coloured with the individual number, ranging from 1 to 420. The black star represents the best individual with the minimum cost function.

Figure 22. Long description

Panel A: A scatter plot shows the evolution of the cost function over evaluations. The x-axis represents the number of evaluations, and the y-axis represents the cost function values. The data points form a downward trend, indicating a decrease in the cost function over time. Panel B: A scatter plot illustrates the lift coefficient increment over evaluations. The x-axis represents the number of evaluations, and the y-axis represents the lift coefficient increment values. The data points show a general upward trend, indicating an increase in the lift coefficient increment over time. Panel C: A scatter plot displays the actuation probability over evaluations. The x-axis represents the number of evaluations, and the y-axis represents the actuation probability values. The data points form a downward trend, indicating a decrease in actuation probability over time. Panel D: A scatter plot shows the individual control laws during LGP. The x-axis represents the change in lift coefficient, and the y-axis represents the actuation probability. The data points are colored according to the individual number, ranging from 1 to 420. The black star represents the best individual with the minimum cost function.

Different from DQN and DDPG, which represent the control policy as a black-box neural network, LGP uses analytical expressions to link the state vector with the control action. By evolving a population of candidate expressions through cross-over and mutation, a satisfying control law can be derived after hundreds of evaluations (Li et al. Reference Li, Noack, Cordier, Borée, Kaiser and Harambat2017b ; Zhou et al. Reference Zhao, Liu, Luo, Deng, Zhou, Peng and Zhang2020; Maceda et al. Reference Maceda, Li, Lusseyran, Morzyński and Noack2021, Reference Maceda, Varon, Lusseyran and Noack2023). For fair comparison, LGP is also implemented in the same hardware as the DRL (CRIO-9038), to see whether better control strategies can be found. The majority of the hyper-parameter settings, such as the generation number and mutation rate, are inherited from Passaggia et al. (Reference Parezanović, Cordier, Spohn, Duriez, Noack, Bonnet, Segond, Abel and Brunton2022). Following Chen et al. (Reference Chen, Wang, Yan, Hu and Noack2025b ), the cost function $J$ is constructed as $J=P_a-\Delta C_l/w_0$ . Namely, two objectives are considered: one to maximise the lift coefficient increment and the other to minimise the mean activation probability (equivalent to actuator power consumption). A weight $w_0$ is divided by the relative lift increment to balance these objectives. To guarantee that the second term has a comparable variation range as the first term ( $P_a$ : 0–1), the weight coefficient $w_0$ is chosen as the maximum value of $\Delta C_l$ (0.148). Other details about the implementation of LGP are proved in Appendix A.

Figure 22(a) plots the optimisation process of the cost function in LGP. For each generation, individuals are ranked in descending order of $J$ . The control laws for the first generations are randomly initialised, resulting in a relatively disperse distribution of $J$ . As the GP proceeds, $J$ decreases steadily and converges after the fifth generation. During the evolutionary process of LGP, the lift increment gradually increases while the actuator power consumption noticeably decreases; see figures 22(b) and 22(c). Figure 22(d) visualises the time evolution of the control laws in a map spanned by $P_a$ and $\Delta C_l$ . Each control law is represented by a dot, with the colour representing its number ( $i$ ). Initially, individuals are widely distributed in the top-left half, with most expressions exhibiting a high energy consumption and low lift increment. After several generations of evolution, the population progressively converges toward the bottom-right corner, finding strategies exhibiting both low energy consumption and high lift gain. The best individual with the minimum cost delivers a relative lift increment of 0.148 and a mean activation probability of 0.145. Correspondingly, the lift coefficient is improved by approximately 21 %. This magnitude of lift improvement is considered as descent for plasma actuators, since its induced jet velocity is very weak, typically between 1–5 ms⁻¹, much less than those of pulsed/synthetic jets.

As summarised in table 3, the lift gain in LGP is identical to the two DRL control algorithms (DDPG and DQN), yet the mean activation probability lies in between, resulting in a slightly lower control efficiency than DDPG (7.71 vs 8.59). In terms of the physical test time, approximately 2 hours are required for LGP to converge, contrasting with the tens of seconds convergence time in DRL. Note that this extremely fast optimisation speed in DRL can be a big deal in practice, since it not only cuts down significantly the operation cost in wind tunnel tests, but also enables real-time tuning of the control law even in flight conditions. Despite the weakness in learning efficiency, the authors did notice that LGP is more robust than DRL, in the sense that DRL requires careful tuning of the reward setting to achieve the best performance. The final control law obtained by LGP reads as

Table 3.

Performance comparison of different closed-loop control methods.

(6.1) \begin{align} \begin{aligned} a_t&=\left [ \cosh \big (\cosh s_3\big ) - \sinh \!\left (\cosh \left (\cos \big ({\rm e}^{{\rm e}^{\sin s_2}}\big )\sinh \left ({\rm e}^{\cosh ^2 s_3}\right )\right )\right ) - \sin \big (\sinh s_4\big ) \right ] \\ &\quad \times \sinh \!\big (\sin \big (\cosh \big (\cosh \big ((s_1 s_4)^2\cos \big (\cosh \big (\cos s_4\big )\big )\big )\big )\big )\big ) .\end{aligned} \end{align}

Figure 23.

A representative state-action time series obtained under (a) DDPG control, (b) DQN control, and (c) LGP control. In each row, the middle plot shows the distribution of action versus state and the right plot shows the probability distribution of the action in different intervals of the state (bin width: 0.6).

To get a taste of how the three closed-loop control strategies behave, the state-action time series obtained by the three strategies are plotted in figure 23. Under DDPG control, a strong phase correlation between the agent’s action and the wake state can be picked. Specifically, the agent is silent when $s_t \lt 0$ , whereas it adopts continuous, high-frequency, wide-amplitude actuation at $s_t \gt 0.6$ . This strategy strives to trap the flow in a high-reward state, which is aligned with the conclusion drawn from the cluster transfer graph (see figure 20). In contrast, under DQN control, the action jumps between 0 and 1, and much more control efforts are invested in the flow separation process ( $s_t \lt 0$ ), resulting in a frequent switch on of the plasma actuation. This ultimately leads to a mean actuation probability approximately three times higher than the DDPG control strategy (0.4 vs 0.13). Under LGP control, the strategy is represented by an analytical expression, and in principle, its output should be continuous. However, practically, a vast majority of the action values are chosen from discrete values of either 0 or 1, and very few fall in between (see the middle plot in figure 23 c). This distribution of the action is very similar to the Bernoulli distribution in the DQN. Different from the latter, plasma actuation in LGP is only shot when the state $s_t$ is located in the vicinity of −1.5 or 2 (see the right plot in figure 23 c), thus yielding a sparse action sequence and a mean actuation probability (0.145) comparable to that of DDPG.

6.2. Turbulent statistics and momentum balance

Figure 24(a) compares the time-averaged velocity profiles extracted at five representative streamwise positions. Under baseline conditions, the main flow is separated from the suction surface at approximately $x=0.13c$ , and at $x=0.15c$ a small portion of reverse flow can already be observed at the bottom of the boundary layer. While moving downstream, the boundary layer velocity profile becomes even more degraded, as indicated by the increasing height of the reverse flow region. The three closed-loop control cases (DDPG, DQN and LGP) introduce almost identical changes to the above flow field, which is largely expected since similar relative lift increments are achieved. To be more specific, near the separation point ( $x=0.2c$ ), closed-loop control barely alters the velocity distribution, indicating that plasma actuation does not energise the inner boundary layer. At $x=0.4c$ , closed-loop control even leads to an increase in the peak reversal velocity at the near-wall region (−0.3 $U_{\infty }$ to −0.5 $U_{\infty }$ ), accelerating the appearance of peak reversal velocity in the recirculation flow. Moving further downstream ( $x\gt 0.6c$ ), the favourable impacts of closed-loop control gradually manifest. Compared with the baseline case, closed-loop control obviously increases the momentum in the bottom half of the boundary layer, resulting in a reduction of the height of the reverse flow. However, this weakened reverse flow in the near-wall region comes at a cost of decreasing velocity in the outer free stream; see figures 24(e) and 24(f). These variations justify that closed-loop plasma actuation suppresses the aerofoil flow separation via momentum redistribution instead of momentum addition.

Figure 24.

(a) Illustration of the five streamwise stations selected to extract the mean velocity profile. The in-plane velocity field under baseline conditions is plotted for reference. (b–f) Comparison of the mean streamwise velocity profiles at $x=0.2c, 0.4c, 0.6c, 0.8c$ and $1c$ .

Figure 24. Long description

Panel A: A contour plot illustrates the in-plane velocity field under baseline conditions. The x-axis is labeled x/c ranging from 0 to 1, and the y-axis is labeled y/c ranging from −0.3 to 0.3. The contour plot shows velocity magnitude with color gradients from light to dark blue, indicating varying velocities. Black lines outline the contour regions, and a black cross marks a specific point of interest. Panel B: A line graph compares the mean streamwise velocity profiles at two different positions. The x-axis is labeled U/U∞ ranging from −0.5 to 1.5, and the y-axis is labeled y/c ranging from −0.3 to 0.3. The graph includes four lines representing different conditions: Baseline (gray), DDPG (blue), DQN (red), and LGP (orange). Panel C: Another line graph shows the mean streamwise velocity profiles at a different position. The axes and lines are labeled similarly to Panel B. Panel D: A line graph depicts the mean streamwise velocity profiles at another position, with similar axes and line representations. Panel E: A line graph illustrates the mean streamwise velocity profiles at yet another position, maintaining consistent axis labels and line representations. Panel F: The final line graph shows the mean streamwise velocity profiles at the last position, with the same axis labels and line representations as the previous panels. Across all panels, the graphs show how velocity profiles change under different conditions and positions, highlighting trends and variations in the flow dynamics.

Figure 25.

Comparison of the second-order turbulent statistics in different cases. In each row, from left to right, the four cases are: (a) baseline, (b) DDPG control, (c) DQN control, and (d) LGP control. In each column, from top to bottom, the quantities displayed are the TKE, streamwise Reynolds normal stress, wall-normal Reynolds normal stress, Reynolds shear stress and turbulence production rate.

For different cases, the spatial distribution of TKE (denoted as $k$ ) is shown in the first row of figure 25, by definition, $k=(\left \langle u'u' \right \rangle +\left \langle v'v' \right \rangle )/2$ . Here $U_{\infty }^2$ is used for normalisation and the resulting dimensionless TKE lies approximately between 0 and 0.2. In the baseline case, a high level of TKE is primarily found in the vicinity of the downstream portion of the shear layer, confined within a streamwise range of $0.5 \lt x/c \lt 1$ . Under closed-loop control, similar spatial redistributions of the TKE are found, irregardless of the difference in control strategy. Specifically, the resultant high-TKE region takes a cone shape that expands downstream from the separation point. The three in-plane Reynolds stresses are further compared in figures 25(e) and 25(f). Overall, the spatial distribution of the streamwise Reynolds normal stress $\left \langle u'u' \right \rangle$ exhibits a strong spatial correlation with the TKE, indicating that the majority of the TKE is caused by the streamwise component. Compared with the baseline case, the most pronounced variation occurs in the transverse Reynolds stress $\left \langle v'v' \right \rangle$ , which spreads widely in the wall-normal direction. Such an increase in the vertical velocity fluctuation is inferred to be caused by the enhanced flapping motion of the separated shear layer under plasma actuation. As with the Reynolds shear stress, the closed-loop control cases exhibits a noticeable increase in both the magnitude and spatial extent of the high shear stress region. According to Pope (Reference Pino, Schena, Rabault and Mendez2000), the turbulence production rate (denoted as $P_k$ ) can be calculated by the relation $P_k=-\langle u'v' \rangle \boldsymbol{\cdot }\partial U/\partial y$ . Consequently, a high-level of turbulence production rate can be found in a narrow band pivoting the centreline of the separated shear layer as a result of the concurrent appearance of the strong velocity gradient and the high-level Reynolds shear stress.

For the purpose of momentum balance analyses, the steady Reynolds-averaged momentum equation along the streamwise direction is written as

(6.2) \begin{equation} U\frac {\partial U}{\partial x}+V\frac {\partial U}{\partial y}+W\frac {\partial U}{\partial z}=-\frac {\partial \langle u^{\prime }u^{\prime }\rangle }{\partial x}-\frac {\partial \langle u^{\prime }v^{\prime }\rangle }{\partial y}-\frac {\partial \langle u^{\prime }w^{\prime }\rangle }{\partial z}+\nu \boldsymbol{\nabla} ^{2}U- \frac {1}{\rho }\frac {\partial p}{\partial x}. \end{equation}

In the present scenario of two-dimensional aerofoil flow separation, the three spanwise gradient terms (i.e. $W \partial U /\partial z$ , $\partial \langle u^{\prime }w^{\prime }\rangle /\partial z$ and $\nu \partial ^{2}U/\partial ^{2}z$ ) can be readily neglected. Moving the cross-stream convection terms to the right-hand side, we can get a simplified version of the streamwise momentum equation:

(6.3) \begin{equation} U\frac {\partial U}{\partial x}=-V\frac {\partial U}{\partial y}-\frac {\partial \langle u^{\prime }u^{\prime }\rangle }{\partial x}-\frac {\partial \langle u^{\prime }v^{\prime }\rangle }{\partial y}+\nu \left (\frac {\partial ^{2}U}{\partial ^{2}y}+\frac {\partial ^{2}U}{\partial ^{2}x}\right )- \frac {1}{\rho }\frac {\partial p}{\partial x} .\end{equation}

Figure 26.

Comparison of the streamwise momentum budget in different cases. In each row, from left to right, the four cases are: (a) baseline, (b) DDPG control, (c) DQN control, and (d) LGP control. In each column, from top to bottom, the displayed quantities are listed on the right-hand side of the colourbar.

Figure 26. Long description

The image contains multiple heatmap graphs comparing the streamwise momentum budget in different cases. Each row represents a different control method: baseline, DDPG control, DQN control, and LGP control. Each column represents different quantities displayed on the right-hand side of the color bar. Panel A: The first row shows the streamwise momentum budget for the baseline case. The x-axis is labeled x/c and the y-axis is labeled y/c. The color bar ranges from −0.5 to 0.5. Panel B: The second row shows the streamwise momentum budget for the DDPG control case. The x-axis is labeled x/c and the y-axis is labeled y/c. The color bar ranges from −0.5 to 0.5. Panel C: The third row shows the streamwise momentum budget for the DQN control case. The x-axis is labeled x/c and the y-axis is labeled y/c. The color bar ranges from 0.5 to 0.5. Panel D: The fourth row shows the streamwise momentum budget for the LGP control case. The x-axis is labeled x/c and the y-axis is labeled y/c. The color bar ranges from −0.5 to 0.5. Each panel shows the distribution of the streamwise momentum budget with red indicating positive values and blue indicating negative values.

Figure 27.

Sketch of the control volume used in momentum budget analysis.

In (6.3), all the velocity-related terms can be derived from PIV measurements. Figure 26 compares the contributions of different terms in (6.3). Note that all the quantities displayed are normalised by $U_\infty ^2/c$ , and in each row the range of the colourbar corresponds to the limits of the normalised quantity. Judging from the magnitudes of different quantities, we can conclude that the two convection terms play a dominant role in the streamwise momentum balance, while the turbulent diffusion and the viscous diffusion terms play secondary roles. In figure 26(ai) the four cases show a similar spatial distribution of the streamwise convection term, characterised by a red positive patch (flow acceleration) located just above the leading edge caused by the suction peak and a blue negative streak (flow deceleration) resulting from the flow separation line. These variations are consistent with the mean velocity profiles shown in figure 24. For the vertical convection term $-V\partial U/\partial y$ , the controlled cases (see figure 26 b) show a noticeable difference when compared with the baseline case. Specifically, the streamwise range of the blue negative patch shrinks obviously, and a large red positive patch takes over the top-right corner. This red positive patch corresponds to the increasing streamwise momentum caused by the downward transportation of the free-stream fluids, serving as direct evidence that plasma actuation can enhance momentum mixing between the two sides of the shear layer. As with the turbulent diffusion terms ( $-\partial \langle u'u'\rangle /\partial x$ and $-\partial \langle u'v'\rangle /\partial y$ ), although their magnitudes are relatively smaller than the mean-flow convection terms, they do play a favourable role in mitigating the flow separation, as evidenced by the increasing level of $-\partial \langle u'u'\rangle /\partial {x}$ just downstream of the separation point (see figure 26 aiii–diii) and the downward shifting of the positive patch of $-\partial \langle u'v'\rangle /\partial y$ (see figure 26 aiv–div). For the viscous diffusion term ( $\nu \boldsymbol{\nabla} ^2 U$ ), no significant difference can be found in the spatial distributions between the baseline and the controlled cases.

To systematically investigate the momentum transport characteristics and the physical mechanisms of flow separation control, a rigorous momentum balance analysis was conducted based on the momentum theorem. As shown in figure 27, the control volume is selected as a quasi-trapezoidal region (denoted as $\varOmega$ ) that fully encloses the suction side of the aerofoil and the entire separated shear layer. The inlet and outlet of the control volume are defined by two vertical lines located at $x = 0c$ and $x = 0.94c$ , respectively. The floor of the control volume is essentially the suction surface of the aerofoil, and the ceiling is closed by a vertical line situated at $y =0.15c$ . For convenience of subsequent notation, the four vertices of the quasi-trapezoidal zone are numbered as A, B, C and D, along the anticlockwise direction. Starting from the conservative form of the momentum equation (6.4) and ignoring the gradient terms along z direction (0 in two-dimensional flow), the following equation can be obtained:

(6.4) \begin{equation} \frac {\partial (\textit{UU})}{\partial x} = -\frac {\partial (\textit{UV})}{\partial y} - \frac {\partial \langle u'u' \rangle }{\partial x} - \frac {\partial \langle u'v' \rangle }{\partial y} + \nu \left ( \frac {\partial ^2 U}{\partial ^2 x} + \frac {\partial ^2 U}{\partial ^2 y} \right ) - \frac {1}{\rho } \frac {\partial p}{\partial x} .\end{equation}

Integrating the first term on the left-hand side within the two-dimensional control volume $\varOmega$ and applying Green’s theorem for simplification results in

(6.5) \begin{align} M_1 = \iint _\varOmega \frac {\partial (\textit{UU})}{\partial x} \mathrm{d}x\mathrm{d}y = \oint _{\partial \varOmega } UU \mathrm{d}y = \left . \int _{y_C}^{y_D} UU \right |_{x=x_D} \mathrm{d}y - \left . \int _{y_B}^{y_A} UU \right |_{x=x_A} \mathrm{d}y ,\end{align}

where $\partial \varOmega$ denotes the boundary of the domain, namely $\overrightarrow {\textit{ABCDA}}$ ; $y_A, y_B, y_C$ and $y_D$ denote the vertical coordinates of the four vertices, and their values are $0.15c, 0c, -0.35c$ and $0.15c$ , respectively.

By analogy, the surface integral of the vertical convection term, turbulent transportation term, viscous diffusion term and pressure gradient term on the right-hand side of (6.4) can be reformulated as

(6.6) \begin{equation} \left \lbrace \begin{array}{l} \displaystyle M_2 = \iint _\varOmega -\frac {\partial (\textit{UV})}{\partial y} \mathrm{d}x\mathrm{d}y = \oint _{\partial \varOmega } UV \mathrm{d}x = -\left . \int _{x_A}^{x_D} UV \right |_{y=y_A} \mathrm{d}x, \\[12pt] \displaystyle M_3 = \iint _\varOmega \left ( -\frac {\partial \langle u'u' \rangle }{\partial x} - \frac {\partial \langle u'v' \rangle }{\partial y} \right ) \mathrm{d}x\mathrm{d}y \\[12pt] \ \ \quad \displaystyle = \left . \int _{y_B}^{y_A} \langle u'u' \rangle \right |_{x=x_A} \mathrm{d}y - \left . \int _{y_C}^{y_D} \langle u'u' \rangle \right |_{x=x_C} \mathrm{d}y - \left . \int _{x_A}^{x_D} \langle u'v' \rangle \right |_{y=y_A} \mathrm{d}x ,\\[12pt] \displaystyle M_4 = \iint _\varOmega \nu \left ( \frac {\partial ^2 U}{\partial ^2 x} + \frac {\partial ^2 U}{\partial ^2 y} \right ) \mathrm{d}x\mathrm{d}y, \\[12pt] \displaystyle M_5 = \iint _\varOmega -\frac {1}{\rho } \frac {\partial p}{\partial x} \mathrm{d}x\mathrm{d}y = -\oint _{\partial \varOmega } \frac {p}{\rho } \mathrm{d}y, \end{array} \right . \end{equation}

where $x_A=0c$ and $x_D=0.94c$ denote the horizontal coordinates of vertices A and D.

Figure 28.

Integral of different terms in the momentum equation for baseline and controlled cases: (a) streamwise convection term $M_1$ , (b) vertical convection term $M_2$ , (c) turbulent transportation term $M_3$ , and (d) viscous diffusion term $M_4$ . Note that the vertical range of the last axis is one percent of the other axes.

Based on the PIV measurement results, we are able to compute almost all the terms in the integral momentum balance equation, except for the pressure gradient term $M_5$ . Figure 28 compares the values of these terms in baseline and controlled cases. For the streamwise convection term, its integral values under all three closed-loop control conditions (DDPG, DQN, LGP) are significantly higher than the baseline case. This can be explained by the fact that plasma actuation effectively suppresses flow separation in the control volume and allows more high-momentum fluids to exit through the right boundary, consistent with the velocity profiles shown in figure 24(f). For the vertical convection term $M_2$ , the baseline case exhibits a large negative integral value, indicating that a significant portion of fluid escapes out of the domain via the ceiling. As a comparison, after plasma actuation, the magnitude of $M_2$ is substantially reduced, demonstrating that closed-loop control indeed enhances the momentum exchange between the shear layer and the main flow, allowing more fluids to reside within the control volume. The turbulent transportation term $M_3$ also plays an important role in the integral momentum balance, and its magnitude varies marginally among different cases, which is largely expected since both baseline and controlled cases exhibit a large-scale unsteady flow separation. The viscous diffusion term $M_4$ plays a marginal role in the momentum balance, and its magnitude ( $O$ (0.0001)) is two orders lower than the other terms. Inter-case comparison shows that the discrete on/off control strategy proposed by DQN is more beneficial than DDPG and LGP in terms of increasing the streamwise momentum within the domain, although their final lift gains are close.

Figure 29.

The distributions of POD modal energy in different cases.

6.3. Proper orthogonal decomposition analyses

Proper orthogonal decomposition has been proven to be an effective approach for identifying the dominant flow structures responsible for the elevated turbulent fluctuations under controlled conditions (Wang et al. Reference Wang, Yan, Hu, Chen, Rabault and Noack2025). For detailed information about the theory, please refer to Berkooz, Holmes & Lumley (Reference Berkooz, Holmes and Lumley1993). Although no dominant frequency information can be extracted, POD analysis remains a popular tool in the community of experimental fluid mechanics since a considerable amount of the experimental data is sparse in the time domain.

Figure 30.

Visualisation of the velocity fields represented by the first four POD modes. From top to bottom, the mode number increases from 1 to 4. From left to right, the four cases displayed are: (a) baseline, (b) DDPG control, (c) DQN control, and (d) LGP control.

Here, the snapshot POD method proposed by Meyer, Pedersen & Özcan (Reference Meyer, Pedersen and Özcan2007) is implemented. The distributions of the modal energy in different cases are plotted in figure 29. Under the baseline condition, the first POD mode alone accounts for approximately 15 % of the total TKE. As the mode number $N_m$ increases, its relative contribution to the total TKE (namely $E_m$ ) drops sharply. At $N_m=2$ , the modal energy is only 7.6 %. Compared with the baseline, the fluctuation energy contained in the first three modes is increased noticeably under closed-loop control. Particularly, the energy contribution of the first POD mode reaches approximately 17.3 % under DRL control and 18.9 % in LGP control. For higher-order POD modes in the range $4\leqslant N_m\leqslant 10$ , their energy shares remain approximately unchanged after plasma actuation.

The velocity fields pertaining to different POD modes are shown in figure 30. According to Zong et al. (Reference Zhu and Krstic2022), the physical meaning of each POD mode can be interpreted by examining the conditional-averaged flow fields constructed from the snapshots that possess the highest or lowest POD mode coefficients. Consequently, the first mode corresponds to the vertical flapping motion of the separated shear layer. Since the area of the positive velocity patch in controlled cases is much larger than that in the baseline case, we can infer that the flapping motion under plasma actuation is more intense, consistent with the elevated modal energy in figure 29. The second mode characterised by a concurrent occurrence of positive velocity patch and negative velocity patch along the streamwise direction, should represent the first-order bending motion of the separated shear layer. Compared with the baseline case, the positive velocity patch is initialised at a more upstream position under controlled conditions, indicating that plasma actuation accelerates the destabilisation process of the separated shear layer. For the third POD mode, a large-scale streamwise vortex can be picked from the spiral velocity pattern above the suction surface. Plasma actuation obviously encourages the formation of this vortex, thereby entraining high-velocity fluids from the free stream into the flow separation zone. The fourth-order mode exhibits different flow patterns in the four cases. Specifically, under DDPG and DQN control, a spiral velocity pattern similar to that in the third POD mode is evident, indicating once again the formation of a streamwise vortex but at a more upstream location. In contrast, the baseline case and the LGP control case seem to correspond to the wall-normal shift of part of the separated shear layer. To sum up, in closed-loop control cases, the high-level TKE beneficial to flow separation suppression is caused by the enhanced flapping motion of the separated shear layer as well as the earlier production of streamwise vortices.