3.1. Convergence history of GCC

Figures 3(a) and 3(b) show the GCC convergence history for the quasiperiodic case. As the number of generations progresses from $j=1$ to $12$ , the number of individuals that achieve low $J$ (below $0.2$ ) increases (figure 3 $a$ ). The $J$ values of the top 50 individuals in each generation are displayed, with the final generation ( $j=12$ ) highlighted in red. The rapid convergence of the GP algorithm is further evident in figure 3( $b$ ), which presents the average $J$ of the top 15 individuals in each generation, denoted as $\bar {J}^j_{1-15}$ . The value of $\bar {J}^j_{1-15}$ quickly settles to a minimum after only the fourth generation. Consequently, for conservativeness, we breed a total of $N_g = 12$ generations for the quasiperiodic case, before extracting the optimal GCC individual, defined as the individual with the lowest $J$ value.

We note that the distribution of $J$ values across the full population need not improve monotonically from one generation to the next. For example, figure 3( $a$ ) shows that generation $7$ contains fewer low- $J$ individuals than generation $5$ , despite the overall convergence trend. This behaviour is expected because, although elitism retains the best-performing individuals, replication, mutation, crossover and the injection of newly generated individuals continue to explore alternative regions of the control-law search space. In this way, intermediate generations may have a broader or less favourable population-level distribution of $J$ , while the best individuals are retained and the evolutionary search continues to converge.

Figure 3.

Convergence history for GCC of quasiperiodic oscillations: ( $a$ ) cost function $J$ of the individuals in each generation and ( $b$ ) the average $J$ of the best 15 individuals in each generation, both shown as the number of generations increases from $j=1$ to $12$ . In ( $c$ ), the control landscape is presented on a 2-D proximity map containing all the GCC individuals, including the optimal one (magenta square) and the uncontrolled case (black triangle). The input–output actuation curve for the optimal GCC individual is depicted in the inset of ( $c$ ).

Figure 3. Long description

Panel A: A line graph shows the cost function of individuals across different generations, with each generation represented by different symbols and colors. The x-axis is labeled ‘Individual’ and ranges from 1 to 50. The y-axis is labeled ‘J’ and ranges from 0 to 1.5. Panel B: Another line graph displays the average cost function of the best 15 individuals in each generation, with the x-axis labeled ‘Generation’ ranging from 1 to 12 and the y-axis labeled ‘J’ ranging from 0 to 1. Panel C: A scatter plot presents the control landscape on a 2-D proximity map, with the x-axis labeled ‘γa’ and the y-axis labeled ‘γb’. The color gradient represents the cost function ‘J’, ranging from 0.16 to 1.4. The optimal individual is marked with a magenta square, and the uncontrolled case is marked with a black triangle. An inset within Panel C shows the input-output actuation curve for the optimal individual.

To visualise the control landscape, we form a 2-D proximity map (figure 3 $c$ ) using classic multidimensional scaling (Young & Householder Reference Young and Householder1938). This gives a centred representation of points $\varGamma =[\gamma _1,\gamma _2,\ldots ,\gamma _{N}]$ with $\gamma _1,\gamma _2,\ldots ,\gamma _{N} \in \mathbb{R}^2$ , such that the pairwise distances between individuals approximate their true distances. To quantify the dissimilarity between two individuals, $K_m$ and $K_n$ (with $1 \leqslant m,n \leqslant N$ ), we define the matrix

(3.1) \begin{equation} D_{m,n}={\big [ \big(B'_{rms,m}-B'_{\textit{rms},n} \big)^2+(\eta _m-\eta _n)^2\big ]}^{1/2}+\alpha |J_m-J_n|, \end{equation}

where $\eta \equiv (p'_{\textit{rms}}-p'_{rms,\ast })/p'_{rms,\ast }$ , with $p'_{\textit{rms}}$ and $p'_{rms,\ast }$ representing the root mean square of the pressure fluctuations with and without control, respectively. The first term in (3.1) captures the difference between two individuals in terms of the actuator power and thermoacoustic amplitude reduction, while the second term penalises disparities in their $J$ values. Following Duriez et al. (Reference Duriez, Brunton and Noack2017), we set $\alpha =2$ to smooth the control landscape. In figure 3( $c$ ), each GCC individual is represented by a discrete marker, whose colour indicates its $J$ value. The uncontrolled case is represented by a black triangle, serving as the reference baseline. The proximity map unveils a scattered pattern, with the best GCC individuals (i.e. those with low $J$ ) congregating towards the far right-hand side, while the remaining individuals are distributed on the left-hand side. The optimal GCC individual exhibits a remarkably low value of the cost function ( $J=0.165$ ), and the functional form of its input–output curve is described by the expression $b(x) = \exp (\sin (1.7x)+\sin (x)-2x)$ . The functional forms of the top 10 GCC individuals and their input–output curves can be found in Appendix A.

3.2. Actuation mechanisms discovered by GCC

Figure 4 shows the system response to the optimal GCC individual from figure 3. Because self-excited thermoacoustic oscillations arise from the two-way coupling between HRR and pressure fluctuations (Rayleigh Reference Rayleigh1878; Lieuwen & Yang Reference Lieuwen and Yang2005), we examine the evolving phase relationship between $q'$ and $p'$ in figure 5. Coupled oscillatory signals can exhibit a range of synchronisation dynamics, including phase drifting, slipping, locking and trapping (Pikovsky et al. Reference Pikovsky, Rosenblum and Kurths2003; Thévenin et al. Reference Thévenin, Romanelli, Vallet, Brunel and Erneux2011). In externally forced self-excited systems, phase trapping and slipping have been reported during transitions from phase drifting to locking (Li & Juniper Reference Li and Juniper2013c ; Balusamy et al. Reference Balusamy, Li, Han, Juniper and Hochgreb2015; Kashinath et al. Reference Kashinath, Li and Juniper2018). In an unforced thermoacoustic system undergoing a quasiperiodic route to chaos, Mondal et al. (Reference Mondal, Pawar and Sujith2017) also identified phase locking, intermittent phase locking and phase drifting between the pressure and HRR signals. These prior observations motivate the analysis below, in which the $p'$ – $q'$ phase relationship is examined in both the uncontrolled state and during the transient and asymptotic stages of GCC. In figures 4 and 5, the evolution of our quasiperiodic system under GCC is presented as five distinct intervals (I–V) according to the stage of control:

Figure 4.

Response of the quasiperiodic system to the optimal GCC individual from figure 3: time traces of the ( $a$ ) acoustic pressure fluctuation $p'$ , ( $c$ ) normalised HRR fluctuation $q'/q'_*$ and ( $e$ ) actuator voltage fluctuation $B'$ . The corresponding PSDs are shown in ( $b$ ), ( $d$ ) and ( $f$ ), respectively.

Figure 5.

Response of the quasiperiodic system to the optimal GCC individual from figure 3: ( $a$ ) evolution of the phase lag between $q'$ and $p'$ , $\Delta \theta _{qp}=\theta _q-\theta _p$ , and ( $b$ ) evolution of the $PLV$ . In interval V, $PLV \approx 1$ , indicating persistence of a synchronised state after control deactivation.

1. (i) Interval I. When uncontrolled (i.e. in its natural state, $0\lt t\lt 2$ s), the system exhibits $p'$ and $q'$ fluctuations at two natural incommensurate frequencies, $f_1$ and $f_2$ , as can be seen in their power spectral densities (PSDs, figures 4 b and 4 d). Notably, the $q'$ mode at the binaural beat frequency ( $f_1-f_2$ ) is particularly strong. In figure 5( $a$ ), the inset ( $0.72\leqslant t \leqslant 0.86$ s) reveals that the phase lag between $q'$ and $p'$ ( $\Delta \theta _{qp}$ ) fluctuates within a bounded range, indicating phase trapping (Aronson et al. Reference Aronson, Ermentrout and Kopell1990). This trapping arises from partial synchronisation between $q'$ and $p'$ , where frequency locking occurs without complete phase locking, which is why this state is sometimes referred to as imperfect phase locking (Li & Juniper Reference Li and Juniper2013c ; Penelet & Biwa Reference Penelet and Biwa2013). We also observe $\Delta \theta _{qp}$ drops in multiples of $2\pi$ radians, known as intermittent synchronisation or phase slipping (Pikovsky et al. Reference Pikovsky, Rosenblum and Kurths2003). Physically, phase slipping signifies that the frequency locking between $q'$ and $p'$ persists for a certain duration until one oscillator momentarily decouples and decelerates, weakening the existing frequency locking (Li & Juniper Reference Li and Juniper2013c ). That oscillator then loses a complete cycle within a short time, manifested as $2\pi$ drops in $\Delta \theta _{qp}$ .

2. (ii) Interval II. Upon activating the optimal GCC individual at $t=2$ s, both $p'$ and $q'$ exhibit time-series modulations resembling beating interactions (figures 4 a and 4 c). During this transient stage, the dominant $p'$ and $q'$ modes maintain their respective natural frequencies (figures 4 b and 4 d). However, a new mode at $f_g=(0.79\pm 0.02) f_1$ appears intermittently in both $p'$ and $q'$ , particularly in the latter. The PSD of the actuation signal ( $B'$ ) is dominated by a strong $f_1$ mode and a weak $f_2$ mode, with little sign of the $f_g$ mode (figure 4 $f$ ). In figure 5( $a$ ), $\Delta \theta _{qp}$ continuously decreases in an unbounded manner. This behaviour, known as phase drifting (Li & Juniper Reference Li and Juniper2013c ), indicates an asynchronous relationship between the two coupled signals. The presence of phase drifting suggests that GCC has disrupted the synchronisation between $q'$ and $p'$ (Mondal et al. Reference Mondal, Pawar and Sujith2017). As observed in the inset (figure 5 $a$ : $3.29\leqslant t \leqslant 3.42$ s), phase drifting in interval II is accompanied by intermittent phase trapping (marked with a dashed square) and occasional phase slipping. The coexistence of phase trapping and drifting signifies that the $q'$ – $p'$ relationship switches temporally between weak synchronisation and desynchronisation.

3. (iii) Interval III. Following the earlier transient stage of beating-like modulations of $p'$ and $q'$ , the amplitude of $p'$ declines gradually, causing the system to settle into a low-amplitude state. In figures 4(b), 4(d) and 4(f), we observe a change in the dominant frequencies of $p'$ , $q'$ and $B'$ , with their natural modes ( $f_1$ and $f_2$ ) giving way to the new mode at $f_g$ . Simultaneously, desynchronisation between $p'$ and $q'$ occurs, as indicated by the continuation of phase drifting (figure 5 $a$ ; inset, $4.42\leqslant t \leqslant 4.53$ s).

4. (iv) Interval IV. The decoupling between $p'$ and $q'$ drives the system towards the final low-amplitude state. Eventually, both $p'$ and $q'$ , along with $B'$ , converge on the $f_g$ mode, although weak signs of the natural modes can still be observed intermittently (figures 4 b, 4 d and 4 f). Notably, the amplitude of $p'$ becomes appreciably lower than its unforced value (figure 4 $a$ ), while the amplitude of $q'$ decreases only slightly (figure 4 $c$ ). In the inset ( $9\leqslant t \leqslant 9.14$ s) of figure 5( $a$ ), we find that $\Delta \theta _{qp}$ remains nearly constant in time, indicating a state of complete synchronisation known as phase locking (Pikovsky et al. Reference Pikovsky, Rosenblum and Kurths2003; Penelet & Biwa Reference Penelet and Biwa2013).

5. (v) Interval V. After GCC is switched off at $t=10$ s, the system remains in the low-amplitude state reached in interval IV (figures 4 a and 4 c). This suggests that GCC has caused the system to evolve to a different attractor (or a different basin of attraction) characterised by low-amplitude $p'$ oscillations. In this post-GCC regime, the residual fluctuations in $p'$ and $q'$ remain strongly phase synchronised, i.e. their phase difference stays nearly constant in time (figure 5 $a$ ). However, this should not be interpreted as the original (natural) thermoacoustic mechanism remaining active. Strong phase synchronisation alone does not imply thermoacoustic growth: self-excited instability requires sufficiently large fluctuation amplitudes and a favourable $p'$ – $q'$ coupling such that the cycle-averaged acoustic power input remains positive. In this post-GCC regime (interval V), the system exhibits coherent but weak fluctuations, consistent with a stabilised low-amplitude synchronised state.

To further examine the phase synchronisation between $q'$ and $p'$ , we compute the phase locking value, $PLV \equiv |\sum _{t=1}^N\exp (i\Delta \theta _{qp})|/N$ , where $\Delta \theta _{qp}$ is the instantaneous phase lag between $q'$ and $p'$ , obtained from the corresponding analytic signals (Lachaux et al. Reference Lachaux, Rodriguez, Martinerie and Varela1999; Mondal et al. Reference Mondal, Pawar and Sujith2017; Roy et al. Reference Roy, Mondal, Pawar and Sujith2020). The phase locking value quantifies the magnitude of the mean complex phase difference and spans the range $[0,1]$ : values near $1$ indicate strong synchronisation with an approximately constant phase difference, whereas values near $0$ indicate weak synchronisation with a widely varying phase difference. Figure 5( $b$ ) shows the phase locking value evolution across intervals I–V. Upon GCC activation at $t=2$ s, the phase locking value decreases and varies strongly, in line with a disruption of the intrinsic coupling and the coexistence of phase trapping and drifting during the transient stage. As the system converges to the controlled low-amplitude state (interval IV), the phase locking value rises towards unity, indicating a constant phase lag between $p'$ and $q'$ . When GCC is deactivated, this phase-locked low-amplitude state persists, with the phase locking value remaining close to unity, consistent with the system remaining on the stabilised attractor reached under GCC.

To understand the origin of the $f_g$ mode, we refer to studies by Roy et al. (Reference Roy, Mondal, Pawar and Sujith2020) and Yin et al. (Reference Yin, Yang, Guan, Redonnet, Gupta and Li2024), who demonstrated that the inherent flame dynamics can play a vital role in determining the frequency at which oscillation quenching occurs in forced self-excited thermoacoustic systems. We therefore conduct additional experiments with the combustor removed, enabling the flame to be examined without acoustic feedback. The open flame is subjected to time-periodic acoustic forcing at a range of frequencies ( $f_f$ ) and velocity-perturbation amplitudes ( $u'/\bar {u}$ ). Here, $u'$ is the velocity fluctuation amplitude at the burner exit, and $\bar {u}$ is the time-averaged bulk velocity. The flame response is characterised by the FDF, which is an amplitude-dependent extension of the flame transfer function (Dowling Reference Dowling1997; Balachandran et al. Reference Balachandran, Ayoola, Kaminski, Dowling and Mastorakos2005; Noiray et al. Reference Noiray, Durox, Schuller and Candel2008). Specifically, the FDF relates the fundamental component of the normalised HRR response at an imposed $f_f$ to the normalised velocity perturbation: $G \equiv (q'/\bar {q})/(u'/\bar {u})$ . Thus, unlike a linear flame transfer function, the FDF permits the gain and phase of the flame response to vary with the forcing amplitude. The FDF approximation assumes that the response at $f_f$ is dominant and that higher harmonics generated by flame nonlinearity have a negligible influence on the overall flame–acoustic dynamics (Dowling Reference Dowling1997; Noiray et al. Reference Noiray, Durox, Schuller and Candel2008).

Figure 6 shows the measured FDF gain. A prominent peak in $G$ emerges at $f_g$ for low values of $u'/\bar {u}$ , indicating that the flame exhibits an exceptionally strong HRR response at this frequency. GCC selects actuation with a spectral component near $f_g$ , exploiting this enhanced receptivity: low-power actuation can then yield a substantial reduction in thermoacoustic amplitude, potentially via AQ (Guan et al. Reference Guan, Gupta, Kashinath and Li2019a ; Mondal et al. Reference Mondal, Pawar and Sujith2019; Roy et al. Reference Roy, Mondal, Pawar and Sujith2020). This finding aligns with our previous work on periodic-regime suppression (Yin et al. Reference Yin, Yang, Guan, Redonnet, Gupta and Li2024), where the GP algorithm autonomously discovered an actuation mechanism that exploits the strong HRR amplification of the open flame near its preferred mode.

For a perfectly premixed flame without equivalence-ratio fluctuations, the velocity-to-HRR transfer-function gain is expected to approach unity in the quasisteady limit $f_f \rightarrow 0$ (Polifke & Lawn Reference Polifke and Lawn2007). Our FDF measurements do not extend to this limit: the lowest forcing frequency is approximately $100$ Hz, which corresponds to $f_f/f_1\approx 0.4$ in the quasiperiodic case. Therefore, the data in figure 6 should not be interpreted as the low-frequency asymptote of the flame response. Nevertheless, within the measured $f_f$ range, GCC identifies an actuation mechanism centred around the preferred mode ( $f_g$ ) of the open flame, despite encountering the flame only in its acoustically coupled configuration.

Figure 6.

Gain of the FDF for an open flame at the quasiperiodic operating condition specified in § 2.1. Experimental data are shown as black dots. The frequency at maximum gain for low $u'/\bar {u}$ is $f_g=(0.79\pm 0.02) f_1$ , which matches the dominant frequency of $B'$ for the optimal GCC individual from figure 3.

The appearance of a spectral component at the preferred (high-gain) mode of the flame under GCC should be interpreted with care. Because the flame exhibits strong receptivity near this mode (figure 6), various actuation waveforms can produce a similar response near this frequency. Thus, the mere presence of $f_g$ components in the controlled spectra does not, by itself, prove that excitation of this mode is the unique causal mechanism of thermoacoustic suppression. Rather, the observation is consistent with GCC converging towards actuation that interacts most effectively with the flame–acoustic system, i.e. near a frequency where the open-FDF has high gain. Crucially, exciting the preferred mode ( $f_g$ ) alone is generally insufficient to suppress thermoacoustic oscillations. The key requirement is to reach and maintain a low-amplitude state with minimal actuation energy. This is where the GP-identified input–output feedback relationship becomes essential. In GCC, the actuation is computed online so that the forcing automatically adapts to the instantaneous system state via pressure measurements. As the thermoacoustic amplitude decreases, the closed-loop mapping typically reduces the actuation effort, enabling the preferred-mode component (when present) to be sustained at low cost while preventing regrowth of the natural thermoacoustic oscillations. This state-dependent adjustment of the actuation signal is consistent with the general operating principle of conventional feedback control, in which the measured system response is used to regulate the actuation demand rather than applying a fixed forcing level (Dowling & Morgans Reference Dowling and Morgans2005). The distinguishing feature of GCC is that its nonlinear feedback map is discovered directly from experimental data by GP, without requiring an explicit model of the thermoacoustic system.

Figure 7.

Instantaneous chemiluminescence snapshots of the flame from the quasiperiodic system at selected time instants: ( $a$ ) uncontrolled state, $t_1$ – $t_5$ , and ( $b$ ) controlled state, $t_6$ – $t_{10}$ . Also shown is ( $c$ ) a comparison of the instantaneous flame fronts between the uncontrolled state (purple) and the controlled state (green), with their corresponding roll-up wrinkle lengths in one representative cycle indicated by $\lambda _*$ and $\lambda$ , respectively.

To investigate the flame dynamics under the influence of the optimal GCC individual, we analyse instantaneous chemiluminescence snapshots at five selected time instants, comparing the uncontrolled and controlled states (figure 7). The flame fronts exhibit characteristic wrinkling caused by travelling waves originating at the flame base (Baillot et al. Reference Baillot, Durox and Prud’Homme1992). Under GCC ( $t_6$ – $t_{10}$ in figure 7 b), the flame fronts oscillate with similar amplitudes to those observed in the uncontrolled state ( $t_1$ – $t_5$ in figure 7 a). This amplitude similarity arises because, while $p'$ is weakened, $q'$ remains comparable in strength between the controlled and uncontrolled states (figure 4 $c$ ). However, the wavelengths of the roll-up wrinkles differ. To illustrate this, we superimpose in figure 7( $c$ ) the flame fronts from the controlled and uncontrolled states. The flame fronts are found by applying an inverse Abel transform to the chemiluminescence images and then extracting the locus of maximum pixel intensity, as per Guan et al. (Reference Guan, He, Murugesan, Li, Liu and Li2019c ). Under GCC, the flame dynamics are governed by the aforementioned preferred mode at frequency $f_g = (0.79 \pm 0.02) f_1$ , which is lower than the dominant frequency $f_1$ of the uncontrolled state. Consequently, the roll-up wrinkles in the controlled state have a longer wavelength ( $\lambda \gt \lambda _*$ ) (Guan et al. Reference Guan, He, Murugesan, Li, Liu and Li2019c ; Roy & Sujith Reference Roy and Sujith2019). Further analysis of the flame dynamics with modal decomposition tools can be found in Appendix C.

3.3. Benchmarking GCC against GOC and OLC

We benchmark the optimal individual from GCC against those from GOC and OLC. For GOC (Duriez et al. Reference Duriez, Brunton and Noack2017; Liu et al. Reference Liu, Tan, Li, Hou, Zhang and Noack2023), the Arduino board used in GCC is replaced with a LabVIEW module (figure 2). The GOC individuals are represented as $b(t)=K[\boldsymbol{x}(t)]$ , where $\boldsymbol{x}(t)$ is a harmonic function capturing the open-loop input dynamics. Specifically, we let $x_i(t)=\sin (2\pi k_if_1t)$ , with $k_i$ values ranging from $0.2$ to $3$ in step sizes of $0.01$ . Figure 8( $a$ ) shows the cost function values ( $J$ ) for six selected GOC generations. The initial generation comprises mostly ineffective individuals, with only a few exhibiting promising control characteristics. Subsequent generations show better performance, as indicated by their lower $J$ values. The average $J$ for each GOC generation rapidly decreases until the third generation, after which it levels off, indicating convergence (figure 8 $b$ ). The optimal GOC individual discovered by the algorithm is

(3.2) \begin{equation} b (t) = \sin [x_1(x_1-x_2)x_3],\,x_i(t)=\sin (2\pi k_if_1t),\,\boldsymbol{k}=[0.61, 0.22, 0.83]. \end{equation}

Figure 8( $c$ ) shows the actuation signal generated by the optimal GOC individual. This signal exhibits repeated cycles, with each cycle being irregular. For OLC, we refer to the results of our previous study (Guan et al. Reference Guan, Gupta, Wan and Li2019b ). The OLC case with the lowest $J$ value is selected as the optimal case for comparison with the optimal GCC and GOC individuals. For brevity, the optimal OLC combination of forcing amplitude and frequency is called the optimal individual as well.

Figure 8.

Convergence history for GOC of quasiperiodic oscillations: ( $a$ ) cost function $J$ of the individuals in each generation and ( $b$ ) the average $J$ of the best eight individuals in each generation, both shown as the number of generations increases from $j=1$ to $10$ . In ( $c$ ), the optimal GOC individual is shown in terms of time traces of the actuation signal.

Figure 9.

Pareto diagrams for the quasiperiodic system under three different control strategies: ( $a$ ) GCC, denoted by filled ‘ $\bigcirc$ ’ markers, ( $b$ ) GOC, denoted by filled ‘ $\vartriangle$ ’ markers, ( $c$ ) OLC, denoted by filled ‘ $\square$ ’ markers and ( $d$ ) all three types of individuals under magnification. The optimal individuals for each control strategy ( $K_{\textit{GCC}}$ , $K_{\textit{GOC}}$ , $K_{\textit{OLC}}$ ) are denoted by hollow magenta markers. Quantitative benchmarking in terms of the cost function $J$ , thermoacoustic amplitude reduction $|\eta |$ , and actuator voltage $B'_{\textit{rms}}$ is shown in ( $e$ ), ( $f$ ) and ( $g$ ), respectively.

Figure 9 uses Pareto diagrams to compare the three control strategies (GCC, GOC, OLC) within a parameter space defined by the input cost ( $B'_{\textit{rms}}$ ) and the state cost ( $\eta$ ). For GCC (figure 9 $a$ ), the most effective individuals are situated in the bottom-left-hand region, achieving a significant reduction in thermoacoustic amplitude while consuming minimal actuation power. Conversely, the least effective individuals are scattered towards the top, with some even amplifying the thermoacoustic amplitude above the uncontrolled reference ( $\eta _{\textit{GCC}} \gt 0$ ). A similar pattern is observed for GOC (figure 9 $b$ ), with the worst individuals concentrated around $\eta _{\textit{GOC}} \approx 0$ , indicating a limited influence on the thermoacoustic amplitude, regardless of the consumed actuation power. As for OLC (figure 9 $c$ ), the best control laws are once again located in the bottom-left-hand region, but the remaining individuals exhibit a wider distribution compared with GCC and GOC, owing to the brute-force nature of open-loop mapping (Guan et al. Reference Guan, Gupta, Wan and Li2019b ). Figure 9( $d$ ) provides a magnified view of the individuals from the three control strategies. The optimal individuals ( $K_{\textit{GCC}}$ , $K_{\textit{GOC}}$ , $K_{\textit{OLC}}$ ), identified by their lowest $J$ values, are highlighted with hollow magenta markers. Although the amplitude reduction achieved by $K_{\textit{GCC}}$ is comparable to that of $K_{\textit{GOC}}$ and $K_{\textit{OLC}}$ , it requires the least actuation effort. Like $K_{\textit{GCC}}$ , both $K_{\textit{GOC}}$ and $K_{\textit{OLC}}$ operate on or near the preferred mode of the open flame (not depicted for brevity) and generate similar thermoacoustic amplitude reductions because of their open-loop nature, which cannot adapt to the specific real-time state of the system. Thus, both $K_{\textit{GOC}}$ and $K_{\textit{OLC}}$ are less efficient than $K_{\textit{GCC}}$ , which incorporates state-dependent feedback. Quantitative benchmarking in terms of $J$ , $|\eta |$ and $B'_{\textit{rms}}$ is shown in figure 9( $e$ – $g$ ). We find that $K_{\textit{GCC}}$ achieves an 87 % reduction in thermoacoustic amplitude with the lowest actuation effort, resulting in the lowest $J$ value of $0.165$ , compared with $J = 0.222$ and $J = 0.260$ for $K_{\textit{GOC}}$ and $K_{\textit{OLC}}$ , respectively.

It is worth noting that the present GCC formulation should not be interpreted as classical linear feedback stabilisation in the narrow sense of continuously stabilising an equilibrium (or a prescribed periodic orbit) via small-signal feedback (Li & Morgans Reference Li and Morgans2016). Instead, the GCC laws provide state-dependent forcing that interacts with the nonlinear flame–acoustic dynamics and can drive the system towards a different dynamical state. In particular, the quasiperiodic example here shows that, once a low-amplitude state is reached, it may persist even after actuation is removed. This observation is consistent with an interpretation in terms of basin-of-attraction (attractor) switching: feedback is essential to reliably and efficiently steer the system into the desired low-amplitude state, while sustained actuation may not always be required to maintain it for the operating conditions considered.

In summary, for the quasiperiodic system, GCC, GOC and OLC can all achieve comparable reductions in thermoacoustic amplitude (around $90\,\%$ ). However, GCC attains this suppression with the lowest actuation cost, resulting in the lowest overall $J$ among all the control strategies tested. Thus, the primary advantage of feedback in this regime is not necessarily a larger thermoacoustic amplitude reduction, but a more favourable performance–economy trade off for efficient control.

4.1. Convergence history of GCC

The GCC convergence history for the suppression of chaotic oscillations is shown in figure 10. Figure 10( $a$ ) presents the sorted $J$ values of the individuals from six selected generations within $j=1$ – $10$ . As the number of generations increases, an increasing number of individuals attain low $J$ values. However, the minimum $J$ value attained by the elite individuals of each generation remains almost unchanged up to the final (10th) generation, marked in red. This shows that GCC discovers the optimal individual relatively quickly, within just the first generation, with subsequent generations uncovering even more low- $J$ individuals. Figure 10( $b$ ) provides additional insights into this convergence trend by showing how $\bar {J}^j_{1-15}$ varies across all 10 generations. The value of $\bar {J}^j_{1-15}$ decreases monotonically until the fifth generation and then fluctuates around a low value thereafter, confirming the rapid convergence of GCC. Considering this and the appearance of new low- $J$ individuals in the last several generations, we execute the GCC algorithm for a total of 10 generations.

Figure 10.

Convergence history for GCC of chaotic oscillations: ( $a$ ) cost function $J$ of the individuals in each generation and ( $b$ ) the average $J$ of the best 15 individuals in each generation, both shown as the number of generations increases from $j=1$ to $10$ . In ( $c$ ), the control landscape is presented on a 2-D proximity map containing all the GCC individuals, including the optimal one (magenta square) and the uncontrolled case (black triangle). The input–output actuation curve for the optimal GCC individual is depicted in the inset of ( $c$ ).

The proximity map, shown in figure 10( $c$ ), reveals a distinct pattern where the low- $J$ individuals tend to cluster tightly towards the left-hand side, with the optimal one (magenta square) situated at the extreme left-hand boundary. Conversely, the high- $J$ individuals are more broadly distributed towards the middle and right-hand side, with the uncontrolled case (black triangle) sitting at the core of the distribution. This suggests that the top-performing individuals share similar actuation characteristics, whereas the bottom-performing individuals are more diverse in their actuation policies, as demonstrated by the latter occupying a larger feature space within the proximity map. The optimal GCC individual is mathematically represented by the functional form $b(x)=\sin (\sin (x)/\exp (x^2)) +\cos (x)/(2.2\exp (x+1.6))+ \cos (\sin (\sin (x)))$ . Detailed information regarding the functional forms of the top 10 GCC individuals, and their input–output curves, can be found in Appendix B.

4.2. Actuation mechanisms discovered by GCC

Figures 11(a), 11(d) and 11(g) present time traces of $p'$ , $q'/q'_{max}$ and $B'$ when the chaotic system is subjected to the optimal GCC individual from figure 10. Figures 11(b), 11(e) and 11(h) show the corresponding PSD evolution, while Figures 11(c), 11(f) and 11(i) show PSD comparisons between the uncontrolled and controlled states. Before GCC activation (uncontrolled), both $p'$ and $q'/q'_{max}$ fluctuate chaotically, as evidenced by their irregular waveforms with intermittent spikes and drops (figures 11 a and 11 d) and broad spectral features with dominant modes at $f_1$ and $f_2$ (figures 11 b and 11 e). However, the $f_2$ mode appears more pronounced than the $f_1$ mode in the HRR spectrum (figure 11 $f$ ), while the opposite trend is observed in the pressure spectrum (figure 11 $c$ ). Upon activating GCC at $t=6$ s, we find a significant reduction in the amplitudes of both $p'$ and $q'/q'_{max}$ (figures 11 a and 11 d), while the amplitude of $B'$ initially jumps and then settles to a lower level (figure 11 $g$ ). Under GCC (controlled state), the corresponding PSDs exhibit weak spectral features, although a clear peak at $f_g=(0.81\pm 0.04) f_1$ can be discerned in the HRR spectra (figures 11 e and 11 f). Upon GCC deactivation at $t=21$ s, both $p'$ and $q'/q'_{max}$ return to their uncontrolled states without hysteresis, resuming their original chaotic dynamics.

Figure 11.

Response of the chaotic system to the optimal GCC individual from figure 10: time traces of the ( $a$ ) acoustic pressure fluctuation $p'$ , ( $d$ ) normalised HRR fluctuation $q'/q'_{max}$ and ( $g$ ) actuator voltage fluctuation $B'$ . The corresponding PSD comparisons are shown in ( $b,c$ ), ( $e,f$ ) and ( $h,i$ ), respectively.

Figure 11. Long description

Panel A: A line graph displays the time trace of acoustic pressure fluctuation (p’) in Pascals over time in seconds. The graph includes a zoomed-in inset highlighting a specific segment of the data. Panel B: A spectrogram shows the power spectral density (PSD) of the acoustic pressure fluctuation, with frequency ratio (f_p/f_1) on the y-axis and time in seconds on the x-axis. The spectrogram highlights different frequency components. Panel C: A line graph compares the PSD of controlled and uncontrolled acoustic pressure fluctuations, with frequency ratio on the y-axis and PSD in dB Hz^-1 on the x-axis. Panel D: A line graph shows the time trace of normalized heat release rate (HRR) fluctuation (q’/q’_max) over time in seconds, including a zoomed-in inset. Panel E: A spectrogram displays the PSD of the normalized HRR fluctuation, with frequency ratio on the y-axis and time in seconds on the x-axis. Panel F: A line graph compares the PSD of controlled and uncontrolled normalized HRR fluctuations, with frequency ratio on the y-axis and PSD in dB Hz^-1 on the x-axis. Panel G: A line graph shows the time trace of actuator voltage fluctuation (B’) in volts over time in seconds, including a zoomed-in inset. Panel H: A spectrogram displays the PSD of the actuator voltage fluctuation, with frequency ratio on the y-axis and time in seconds on the x-axis. Panel I: A line graph compares the PSD of controlled actuator voltage fluctuations, with frequency ratio on the y-axis and PSD in dB Hz^-1 on the x-axis.

Figure 12 presents the synchronous behaviour between $q'$ and $p'$ . The entire process is divided into four intervals, labelled as I–IV. In interval I (uncontrolled state), when GCC is not yet applied, $\Delta \theta _{qp}$ decreases extremely gradually, with short periods of constant phase lag appearing intermittently (figure 12; inset, $2.1\leqslant t \leqslant 2.4$ s). This is consistent with intermittent phase drifting and suggests desynchronisation, accompanied by brief episodes of phase locking. Upon GCC activation at $t=6$ s, $q'$ and $p'$ become phase locked, as evidenced by the constant phase lag observed in interval II. However, after a short stint, $\Delta \theta _{qp}$ begins to rise, leading to phase drifting in interval III. The GCC forcing disrupts the coupling between $p'$ and $q'$ , facilitating their desynchronisation, until GCC deactivation at $t=21$ s. Subsequently, the system returns to its original uncontrolled state of chaos in interval IV.

Figure 12.

Response of the chaotic system to the optimal GCC individual from figure 10: evolution of the phase lag between $q'$ and $p'$ , $\Delta \theta _{qp}=\theta _q-\theta _p$ .

In a manner similar to the quasiperiodic case, we examine the FDF at the chaotic operating condition in figure 13. As described in § 3.2, the FDF characterises the fundamental HRR response of an open flame to harmonic velocity perturbations as a function of the forcing frequency and amplitude. It is measured with the combustor removed and thus represents the flame-only response rather than the coupled thermoacoustic dynamics. We repeat the FDF measurements here because the chaotic case has a different bulk reactant velocity and equivalence ratio from the quasiperiodic case (§ 2.1).

For low forcing amplitudes, figure 13 shows a prominent gain peak at $f_g$ , corresponding to the preferred mode of the open flame (Roy et al. Reference Roy, Mondal, Pawar and Sujith2020; Yin et al. Reference Yin, Yang, Guan, Redonnet, Gupta and Li2024). The elevated value of $G$ indicates that $q'$ is highly receptive to perturbations near this frequency. The preferred mode coincides with the weak spectral component at $f_g$ observed in the GCC state (figure 11). As in the quasiperiodic case, the FDF measurements begin at a finite forcing frequency of approximately $100$ Hz and therefore do not resolve the quasisteady limit in which the gain of a perfectly premixed flame without equivalence-ratio fluctuations is expected to approach unity (Polifke & Lawn Reference Polifke and Lawn2007).

The presence of a preferred mode in the controlled spectra is consistent with GCC exploiting a high-receptivity spectral pathway, but does not by itself establish a unique causal mechanism. The decisive element is the state-dependent feedback actuation provided by GCC, which adapts the forcing to the chaotic dynamics while penalising the actuation effort. Thus, the GCC algorithm effectively combines the enhanced HRR response near the preferred mode with state-dependent feedback control, converting low-power actuation into a marked reduction in thermoacoustic amplitude.

Figure 13.

Gain of the FDF for an open flame at the chaotic operating condition specified in § 2.1. Experimental data are shown as black dots. The frequency at maximum gain for low $u'/\bar {u}$ is $f_g=(0.81\pm 0.04) f_1$ , which matches the preferred-mode frequency from figure 11.

The influence of GCC is visually examined in figure 14 via instantaneous flame chemiluminescence snapshots at selected times, for both the uncontrolled and controlled states. In the uncontrolled state (figure 14 a, $t_1$ – $t_5$ ), the flame fronts exhibit pronounced surface wrinkling caused by a travelling wave running along the flame body (Baillot et al. Reference Baillot, Durox and Prud’Homme1992). The lift-off height of these chaotic flames varies across different time intervals. In the controlled state (figure 14 b, $t_6$ – $t_{10}$ ), the flame fronts appear smoother with only minor low-amplitude movement. This enhanced steadiness is a result of the amplitudes of both $p'$ and $q'$ being broadly reduced to a low level (figures 11 a and 11 d). Further analysis of the flame dynamics can be found in Appendix D.

Figure 14.

Instantaneous chemiluminescence snapshots of the flame from the chaotic system at selected time instants: ( $a$ ) uncontrolled state, $t_1$ – $t_5$ , and ( $b$ ) controlled state, $t_6$ – $t_{10}$ .

4.3. Characterising the controlled state under GCC

In § 4.2, we demonstrated that under GCC, the chaotic system exhibits an exceedingly low $p'$ amplitude, with its PSD displaying weak broadband characteristics. We now examine these $p'$ dynamics in more detail, focusing on distinguishing between chaotic and stochastic dynamics in the controlled state.

We first invoke the 0–1 test from Gottwald & Melbourne (Reference Gottwald and Melbourne2004). This test can distinguish between regular and chaotic dynamics in deterministic systems (Guan et al. Reference Guan, Yin, Yang and Li2024). Chaotic dynamics yield a test value near unity, whereas regular dynamics yield a test value near zero. Figure 15( $a$ ) shows the motion of the translation variables ( $m_c,n_c$ ) from the 0–1 test, which are cumulative sums of the input time series after projecting it onto two reference patterns. A visual examination reveals that both the uncontrolled and controlled states host Brownian-like trajectories, indicating chaos. Although not shown for brevity, the 0–1 test values for both states are approximately unity, confirming the presence of chaos in both states.

Figure 15.

Time-series analysis of the $p'$ signal for the uncontrolled and controlled states of the chaotic system. Panel ( $a$ ) shows the translation variables ( $m_c,n_c$ ) from the 0–1 test, while ( $b$ ) shows the BP probability distributions, all using $20\,000$ -point time-series data (each $2.0$ s in duration) and an embedding dimension of $d=5$ . Panel ( $c$ ) shows the CECP: permutation statistical complexity ( $C_{JS}$ ) versus permutation entropy ( $H$ ), computed from $5000$ -point time-series data ( $0.5$ s in duration), with $d=6$ . In ( $c$ ), the black solid curves denote the theoretical upper and lower bounds of $C_{JS}$ , and the blue dashed curve corresponds to fractional Brownian motion (fBm) (Mandelbrot & Van Ness Reference Mandelbrot and Van Ness1968), included here to assess stochasticity.

To corroborate this, we examine the Bandt–Pompe (BP) probability distribution based on the permutation spectrum (Bandt & Pompe Reference Bandt and Pompe2002; Miyano & Gotoda Reference Miyano and Gotoda2017; Zhu et al. Reference Zhu, White, Carter, Baek and Terry2017). This spectrum quantifies the variability of ordinal pattern spectra across different disjoint segments of an input time series (Amigó et al. Reference Amigó, Kocarev and Szczepanski2006; Kulp & Zunino Reference Kulp and Zunino2014; Guan et al. Reference Guan, Murugesan and Li2018b ). We use an embedding dimension of $d=5$ to reconstruct permutations of the ordering of $p'$ amplitudes from consecutive time-series elements. This produces ( $N-d+1$ ) $d$ -tuples ( $a_1, a_2, \ldots , a_{N-d+1}$ ) from a time series with $N$ points. The amplitude ordering combinations have $d!$ possible permutations (Riedl et al. Reference Riedl, Müller and Wessel2013; Tony et al. Reference Tony, Gopalakrishnan, Sreelekha and Sujith2015). The BP probability measures the occurrence probability of structures within the time series. If $ \pi _1, \ldots , \pi _{d!}$ denotes all possible permutation patterns, the occurrence probability of each pattern can be calculated as $P_m=n(\pi _m)/(N-d+1)$ (Zhu et al. Reference Zhu, White, Carter, Baek and Terry2017). A chaotic signal contains structures with preferred states in the probability space, whereas a stochastic signal exhibits a relatively uniform BP distribution, lacking a preferred amplitude ordering and occupying the entire BP probability space with equal likelihood (Bandt & Pompe Reference Bandt and Pompe2002). In our GCC experiments (figure 15 $b$ ), both the uncontrolled and controlled states exhibit non-uniform BP probability distributions, with elevated values at low bin numbers, indicating some degree of chaos in both states.

The shape of the BP probability distribution offers qualitative insights into the nature of a given time series, and this can be further quantified using the complexity–entropy causality plane (CECP) (Rosso et al. Reference Rosso, Larrondo, Martin, Plastino and Fuentes2007a , Reference Rosso, Zunino, Pérez, Figliola, Larrondo, Garavaglia, Martín and Plastinob ). The CECP is a 2-D map formed by plotting the permutation statistical complexity ( $C_{JS}$ ) against the permutation entropy ( $H$ ). Based on the BP probability ( $P_m$ ), for a given probability distribution $\boldsymbol{P}=(P_1, P_2, \ldots , P_{d!})$ , $H$ can be expressed in terms of the normalised Shannon entropy $S$ (Bandt & Pompe Reference Bandt and Pompe2002):

(4.1) \begin{equation} H(\boldsymbol{P})=\frac {S(\boldsymbol{P})}{S_{max}}=-\frac {\sum _{m=1}^{d!}P_m\log _2 P_m}{\log _2d!}. \end{equation}

The statistical complexity of a given probability distribution can be expressed as

(4.2) \begin{equation} C_{JS}(\boldsymbol{P})=Q_{JS}(\boldsymbol{P})\boldsymbol{\cdot }H(\boldsymbol{P}), \end{equation}

where $Q_{JS}(\boldsymbol{P})$ is the disequilibrium found from the Jensen–Shannon divergence ( $D_{JS}$ ) and a normalisation constant $Q_0=-2/[(({d!+1})/{d!})\log _2(d!+1)-2\log _2(2d!)+\log _2(d!)]$ . Thus, $Q_{JS}(\boldsymbol{P})$ can be expressed as

(4.3) \begin{equation} Q_{JS}(\boldsymbol{P})=Q_0\boldsymbol{\cdot }D_{JS}=Q_0\left [S \left(\frac {\boldsymbol{P}+\boldsymbol{P_u}}{2} \right)-\frac {S(\boldsymbol{P})}{2}-\frac {S(\boldsymbol{P_u})}{2}\right ]\!, \end{equation}

where $\boldsymbol{P_u}$ is a uniform distribution. Different regions on the CECP correspond to different types of dynamics. Periodic systems exhibit low entropies and low complexities owing to their repetitive patterns and strong predictability. Chaotic systems, by contrast, exhibit moderate entropies and high complexities. Stochastic systems exhibit high entropies and low complexities owing to their high uncertainties and uniform probability distributions (Soriano et al. Reference Soriano, Zunino, Larger, Fischer and Mirasso2011; Zunino et al. Reference Zunino, Soriano and Rosso2012; Zhu et al. Reference Zhu, White, Carter, Baek and Terry2017). In the CECP, points on or below the fBm line can be considered purely stochastic (Mandelbrot & Van Ness Reference Mandelbrot and Van Ness1968). In our GCC experiments, both the uncontrolled and controlled states of the initially chaotic system are analysed with multiple time delays and plotted on the CECP (figure 15 $c$ ). The black solid curves denote the theoretical upper and lower bounds of $C_{JS}$ , and the blue dashed curve corresponds to fBm. In figure 15( $c$ ), we find that a large number of data points (red markers), representing the uncontrolled state, reside near the upper-left-hand region of the CECP, indicating strong chaos (Weygand & Kivelson Reference Weygand and Kivelson2019; Lee Reference Lee2022). Conversely, the data points for the controlled state (green markers) are positioned towards the right-hand side of the CECP, indicating chaos with strong noise (Soriano et al. Reference Soriano, Zunino, Larger, Fischer and Mirasso2011; Zunino et al. Reference Zunino, Soriano and Rosso2012; Gotoda et al. Reference Gotoda, Kobayashi and Hayashi2017; Murayama et al. Reference Murayama, Kinugawa, Tokuda and Gotoda2018; Hachijo et al. Reference Hachijo, Masuda, Kurosaka and Gotoda2019). This distinction is consistent with our experimental findings, where the $p'$ amplitude is extremely low under GCC (figure 11 $a$ ), making it more susceptible to noise. In summary, based on the 0–1 test and complexity–entropy diagnostics, we conclude that even under GCC-induced thermoacoustic suppression, the chaotic system retains some degree of deterministic chaos.

4.4. Benchmarking GCC against GOC and OLC

The GOC convergence history for the chaotic case is shown in figure 16. As the number of generations progresses from $j=1$ to $10$ , the control performance, as quantified by $J$ , remains largely unchanged (figure 16 $a$ ). This is confirmed by $\bar {J}^j_{1-15}$ (figure 16 $b$ ), which fluctuates around a moderate value (around $0.6$ – $0.8$ ) with increasing $j$ . Relative to GCC, the inferiority of GOC is attributed to the absence of state-dependent feedback that continually adjusts the actuation in real time based on the system dynamics. The GOC actuation command is predetermined and applied without considering the actual response of the system, which can lead to suboptimal performance and even amplified oscillations in some cases.

Figure 16.

Convergence history for GOC of chaotic oscillations: ( $a$ ) cost function $J$ of the individuals in each generation and ( $b$ ) the average $J$ of the best 15 individuals in each generation, both shown as the number of generations increases from $j=1$ to $10$ .

To quantitatively evaluate the performance of the optimal GCC individual in comparison with GOC and OLC, we consider the cost function $J$ , thermoacoustic amplitude reduction $\eta$ and actuation voltage $B'_{\textit{rms}}$ , as shown in figure 17. For GCC, most of the individuals congregate towards the left-hand side of the Pareto diagram (figure 17 $a$ ), where the required actuation effort is low. The top-performing GCC individuals are concentrated at the bottom of this left-hand band, where the amplitude reduction is high. The optimal GCC individual ( $K_{\textit{GCC}}$ ), denoted by a magenta ‘ $\bigcirc$ ’ marker in figure 17( $a$ ), achieves the highest amplitude reduction ( $|\eta | = 95.7\,\%$ ) with the lowest actuation effort ( $B'_{\textit{rms}} = 0.03$ V), yielding the lowest cost function ( $J = 0.082$ ), as shown in the benchmarking plots of figure 17( $e$ – $g$ ).

Figure 17.

Pareto diagrams for the chaotic system under three different control strategies: ( $a$ ) GCC, denoted by filled ‘ $\bigcirc$ ’ markers, ( $b$ ) GOC, denoted by filled ‘ $\vartriangle$ ’ markers, ( $c$ ) OLC, denoted by filled ‘ $\square$ ’ markers and ( $d$ ) all three types of individuals under magnification. The optimal individuals for each control strategy ( $K_{\textit{GCC}}$ , $K_{\textit{GOC}}$ , $K_{\textit{OLC}}$ ) are denoted by hollow magenta markers. Quantitative benchmarking in terms of the cost function $J$ , thermoacoustic amplitude reduction $|\eta |$ , and actuator voltage $B'_{\textit{rms}}$ is shown in ( $e$ ), ( $f$ ) and ( $g$ ), respectively.

For GOC, the individuals are more scattered across the Pareto diagram (figure 17 $b$ ), indicating greater variability in actuation pathway and control performance. A considerable number of GOC individuals produce $\eta _{\textit{GOC}} \gt 100\,\%$ , implying amplification of the thermoacoustic oscillations relative to the uncontrolled case. The optimal GOC individual ( $K_{\textit{GOC}}$ ), denoted by a magenta ‘ $\vartriangle$ ’ marker in figure 17( $b$ ), cannot match the amplitude reduction of $K_{\textit{GCC}}$ despite using more actuation effort ( $B'_{\textit{rms}} = 0.04$ V), resulting in a higher cost function ( $J = 0.274$ ). Here, $K_{\textit{GOC}}$ does not target the preferred mode of the flame.

For OLC, the individuals are predetermined by brute-force open-loop mapping and occupy a C-shaped region in the Pareto diagram (figure 17 $c$ ). The optimal OLC individual ( $K_{\textit{OLC}}$ ), denoted by a magenta ‘ $\square$ ’ marker in figure 17( $c$ ), can achieve a similar amplitude reduction to $K_{\textit{GCC}}$ , but requires significantly more actuation effort ( $B'_{\textit{rms}} = 0.07$ V), resulting in a near doubling of the cost function ( $J = 0.172$ ) relative to GCC.

All three optimal individuals ( $K_{\textit{GCC}}$ , $K_{\textit{GOC}}$ , $K_{\textit{OLC}}$ ) are superimposed in figure 17( $d$ ). The superiority of $K_{\textit{GCC}}$ over both $K_{\textit{GOC}}$ and $K_{\textit{OLC}}$ is clear: $K_{\textit{GCC}}$ not only produces the largest thermoacoustic amplitude reduction but also uses the least actuation power.

In summary, when applied to the chaotic system, GCC provides the best performance across all evaluation metrics: it achieves the largest reduction in thermoacoustic amplitude while simultaneously requiring the least actuation effort. In the chaotic regime, GCC not only reduces the actuation power needed for a fixed suppression level, but it also improves suppression and efficiency concurrently.