2.1. Stochastic self-sustained oscillator model

Here, we describe briefly our system model and the corresponding probabilistic solution. One may refer to Lee et al. (Reference Lee, Kim and Park2023b) for a more detailed mathematical description. First, we introduce a phenomenological low-order model for a thermoacoustically oscillating system. Specifically, we consider an SDE consisting of a Van der Pol type self-sustained oscillator equation and an additive noise term:

(2.1)\begin{equation} \frac{\mathrm{d}^2x}{\mathrm{d}t^2} - (\epsilon+\alpha x^2)\,\frac{\mathrm{d}\kern0.06em x}{\mathrm{d}t} + \omega^2 x=\sqrt{2d}\,\eta,\quad \text{for }t > 0, \end{equation}

where $x$ is the system variable (e.g. pressure in a combustor), $\eta$ is a unit white Gaussian noise, and $d>0$ is the amplitude of the noise. Here, $\epsilon$ and $\alpha$ are linear and nonlinear parameters, equivalent to $k_1$ and $k_2$ in (1.1), respectively, with scale factors. Finally, $\omega$ is the angular frequency of the oscillation that can be obtained easily from spectral analysis in practice. A probabilistic solution of (2.1) in the form of the Fokker–Planck equation can be obtained by applying the method of variation of parameters. Specifically, we transform the system variable $x$ and its derivative $\mathrm {d}\kern0.06em x/\mathrm {d}t$ into the functions of amplitude ($a$) and phase ($\phi$):

(2.2)\begin{equation} \left.\begin{gathered} x = a\cos{(\omega t+\phi)},\\ \frac{\mathrm{d}\kern0.06em x}{\mathrm{d}t} =-a\omega \sin{(\omega t+\phi)}. \end{gathered}\right\} \end{equation}

The method of variation of parameters shown above is used frequently for the analysis of stochastic nonlinear oscillators (Nayfeh Reference Nayfeh1981; Zhu & Yu Reference Zhu and Yu1987). By applying trigonometric identities, a set of SDEs is obtained:

(2.3)\begin{equation} \left.\begin{gathered} \frac{\mathrm{d}a}{\mathrm{d}t} = \frac{\epsilon}{2}\,a + \frac{\alpha}{8}\,a^3 + Q_1 - \left( \frac{\sqrt{2d}}{\omega} \sin \varPhi \right) \eta_1, \\ \frac{\mathrm{d}\phi}{\mathrm{d}t} = Q_2 - \left( \frac{\sqrt{2d}}{\omega a} \cos \varPhi \right) \eta_2, \end{gathered}\right\} \end{equation}

where $\varPhi$ is $\omega t + \phi (t)$, and $Q_1$ and $Q_2$ are the sums of first-order terms that become zero upon time averaging. Here, $\eta _1$ and $\eta _2$ are unit white Gaussian noise terms for amplitude and phase, respectively, each of which is an independent stochastic process with a correlation time smaller than the acoustic period (Bonciolini et al. Reference Bonciolini, Faure-Beaulieu, Bourquard and Noiray2021; Indlekofer et al. Reference Indlekofer, Faure-Beaulieu, Dawson and Noiray2022). Therefore, we can confirm the equivalence of (2.1), (2.3) and (1.1) for zero noise ($d=0$) when averaged over sufficient time. Applying the stochastic averaging under the assumption of weak nonlinearity, we obtain the Fokker–Planck equation:

(2.4)\begin{align} \left.\begin{gathered} \partial_t P(a,t) =-\partial_a(D^{(1)}(a)\,P(a,t))+\partial_{aa}(D^{(2)}P(a,t)),\quad \text{for } (a,t)\in (0,\infty]\times[0,\infty],\\ P(a,0) = P_0(a),\quad \text{for } a\in[0,\infty], \\ P(0,t) = 0,\quad \text{for } t\in[0,\infty], \end{gathered}\right\}\end{align}

where $D^{(1)}(a) = ({\epsilon }/{2})a + ({\alpha }/{8})a^3 + {d}/{2\omega ^2 a}$ and $D^{(2)}(a) = {d}/{2\omega ^2}$.

2.2. A PINN for solving the inverse problem

In the system identification of thermoacoustic oscillators, we aim to determine the unknown system parameters $\epsilon$, $\alpha$ and $d$ in the SDE (2.1) from the discrete stochastic samples $\{x^{(i)}_{t_j}\}_{i=1}^N$ at various times $t_j$ for $j=1, 2, \ldots, M$. Although recent studies (Chen et al. Reference Chen, Yang, Duan and Karniadakis2021; Xu & Darve Reference Xu and Darve2021; Shin & Choi Reference Shin and Choi2023) have demonstrated successfully data-driven methods for solving such SDE-related inverse problems, training neural networks while adhering to SDE constraints typically demands a considerable number of samples, and may lead to unstable training. In our approach, we tackle these challenges by adopting deterministic surrogate modelling, specifically by integrating a smooth initial condition and domain truncation derived from the Fokker–Planck equation (Son & Lee Reference Son and Lee2023). The surrogate equation for a PINN-based system identification is

(2.5)\begin{equation} \left.\begin{gathered} \partial_t P(a,t) =-\partial_a(D^{(1)}(a)\,P(a,t))+D^{(2)}\partial_{aa}P(a,t), \quad\text{for } (a,t)\in [0,A]\times[0,T], \\ P(a,0) = \hat{P}(a), \quad\text{for } a\in[0,A], \\ P(0,t) = 0,\quad \text{for } t\in[0,T], \end{gathered}\right\}\end{equation}

where $D^{(1)}(a)$ and $D^{(2)}(a)$ are as in § 2.1, and $\hat {P}(a)$ is a gamma distribution $\varGamma (\theta _1, \theta _2)$ with $\theta _1 \geq 2$ and $\theta _2 \gg 1$.

We aim to approximate simultaneously the solution of (2.5) and the unknown model parameters based on (2.5) with the stochastic samples $\{x^{(i)}_{t_j}\}_{i=1}^N$. To achieve this, we utilize a fully connected neural network to approximate the solution, and three additional learnable parameters $\epsilon$, $\alpha$ and $d$ to approximate the unknown model parameters. The network comprises three hidden layers, each composed of 64 hidden units activated by the hyperbolic tangent function. The output layer is fed to the softplus function $F(x) = \log (1+\textrm {e}^x)$, ensuring that the network output remains non-negative. We define three loss functions: $\mathcal {L}_R$ for the PDE residual, $\mathcal {L}_{BC}$ for the boundary condition, and $\mathcal {L}_{mass}$ to ensure that the solution conforms to the properties of a probability density function (p.d.f.). Mathematically, these functions are represented as

(2.6)\begin{equation} \left.\begin{gathered} \mathcal{L}_R = \Vert\partial_t P(a,t) + \partial_a(D^{(1)}(a)\,P(a,t)) - D^{(2)}\partial_{aa}P(a,t) \Vert_{L^2([0,A]\times[0,T])}^2,\\ \mathcal{L}_{BC} = \Vert P(a,t)\Vert_{L^2(\{0\}\times[0,T])}^2, \\ \mathcal{L}_{mass} =\left\Vert \int_{[0,A]} P(a,t)\,{\rm d}a - 1 \right\Vert_{L^2([0,T])}^2, \end{gathered}\right\} \end{equation}

where the learnable parameters $\epsilon$, $\alpha$ and $d$ contribute to $\mathcal {L}_R$ through $D^{(1)}(a)$ and $D^{(2)}$. Note that we excluded the surrogate initial condition from the loss function, as the initial data will be given repeatedly in stochastic samples.

System identification of (2.5) requires integrating the stochastic samples $\{x^{(i)}_{t_j}\}_{i=1}^N$ with the solution of the Fokker–Planck equation (2.5), which represents the p.d.f. of the amplitude variable $a$ at time $t$. Our process begins by transforming the samples $\{x^{(i)}_{t_j}\}_{i=1}^N$ into amplitudes $\{a^{(i)}_{t_j}\}_{i=1}^N$ through the use of the Hilbert transform. To facilitate the integration, we make an independence assumption about the samples $\{x_{t_j}^{(i)}\}_{i=1}^N$ for each time $t_j$ for $j=1,\ldots,M$. This assumption allows us to incorporate the concept of maximum likelihood estimation into the loss function of PINN. Because maximizing likelihood is equivalent to minimizing negative log-likelihood, we define a data loss function by the sum of negative log-likelihood as

(2.7)\begin{equation} \mathcal{L}_{data} = \sum_{j=1}^{N} -\log \left(\prod_{i=1}^{N}P(a_i, t_j)\right) =-\sum_{i,j=1}^{N,M}\log P(a_i, t_j). \end{equation}

Finally, we solve an optimization problem:

(2.8)\begin{equation} \min_{W, \epsilon, \alpha, d} \lambda_1\mathcal{L}_{R} + \lambda_2\mathcal{L}_{BC} + \lambda_3\mathcal{L}_{mass} + \lambda_4\mathcal{L}_{data}, \end{equation}

for some $\lambda _1, \lambda _2, \lambda _3, \lambda _4$, where $W$ denotes the network parameters. Adam optimizer (Kingma & Ba Reference Kingma and Ba2014) is used for the optimization. The overall architecture of the proposed PINN-based system identification is illustrated in figure 1.

Figure 1.

The architecture of the PINN designed for identifying the parameters of the self-sustained thermoacoustic oscillator. We utilize a fully connected neural network providing non-negative output value, and three additional learnable parameters $\epsilon$, $\alpha$ and $d$. Network output is processed through automatic differentiation and combined with $\epsilon$, $\alpha$ and $d$ to compute the loss functions as defined in (2.6) and (2.7). We update both the network parameters and the additional learnable parameters through back propagation until convergence is achieved.

3.1. Synthetic data

For the validation of the PINN-based system identification framework, we use both the synthetic data and the experimental data. First, we consider six sets of synthetic data, each containing ten numerically generated time series. These time series feature stochastic transient self-sustained oscillation displayed in figures 2(a–f). Three sets of data are in the fixed-point regime before the Hopf bifurcation ($\epsilon = -0.3, -0.2, -0.1$), while the other three sets are in the limit-cycle regime after the Hopf point ($\epsilon = +0.1, +0.2, +0.3$). Other system parameters are are set as $\alpha =-0.1$, $d=0.1$ and $\omega =2{\rm \pi}$, as per Son & Lee (Reference Son and Lee2023).

Figure 2.

(a–f) Sample time series (black line) and the amplitude (blue line) used for PINN training. (g–l) Time evolution of the amplitude in ten sample time series at each case, and (m–r) their p.d.f.s. Six sets of synthetic data are displayed: (a,g,m) $\epsilon =-0.3$, (b,h,n) $\epsilon =-0.2$, (c,i,o) $\epsilon =-0.1$ at the fixed-point regime, and (d,j,p) $\epsilon =+0.1$, (e,k,q) $\epsilon =+0.2$, (f,l,r) $\epsilon =+0.3$ at the limit-cycle regime. Here, $\alpha$, $d$ and $\omega$ are fixed at $-0.1$, 0.1 and $2{\rm \pi}$ in all cases.

Synthetic data are created by solving (2.1) numerically with the fourth-order Runge–Kutta method for $t=[0,100]$ with $\mathrm {d}t=0.01$. When solving the time-marching problem, We set the initial values of [$x, \mathrm {d}\kern0.06em x/\mathrm {d}t$] as $[1,0]$ for fixed-point data, and $[0,0]$ for limit-cycle data. These initial conditions enable the observation of transient amplitude death in the fixed-point regime and amplitude growth in the limit-cycle regime, as depicted in figures 2(a–l). Time evolution of the oscillation amplitudes computed from the Hilbert transform are shown collectively in figures 2(g–l), while their p.d.f.s are shown in figures 2(m–r). The saturation in deterministic oscillation amplitude is found at $t>80$ for all synthetic data, indicating zero-amplitude fixed points with stochastic fluctuation ($\epsilon <0$) and fully developed self-sustained oscillations ($\epsilon >0$). It is worth noting that the amplitude saturation is slower in the weakly nonlinear time series where the absolute value of the linear growth rate is close to zero.

3.2. Experimental data

As for the experimental validation, we use the pressure data obtained from an annular model gas-turbine combustor identical to that in Guk et al. (Reference Guk, Seo and Lee2023) (figure 3a). In this set-up, gaseous methane (purity 95.95 %) is premixed with air, which is passed through the dryer (Kyungwon T15K) and compressor (Kyungwon AL5N), before entering the combustor via a swirler and a nozzle. The swirler and the nozzle, respectively, have swirl number 0.608 and diameter 35 mm. The annular combustor has inner diameter 395 mm and outer diameter 405 mm, and is manufactured with SUS304 stainless steel except for the inner wall. A cylindrical quartz serves as the inner wall of the annular combustor, enabling the visual inspection of the flame via two planar mirrors. Nine circular openings with diameter 60 mm are installed at the combustor ceiling, which constitutes the flow exit along with a fan-shaped opening at the top of the nozzle. The methane–air equivalence ratio is varied between 0.74 (lean blowoff limit) and 0.9, with step size 0.02, while the total mass flow rate is kept to $5.5\ \textrm {g}\ \textrm {s}^{-1}$ via mass flow controllers (MKP TSC-230 and MKP TSC-145 for fuel and air, respectively). Pressure oscillation in the combustor is measured with a piezoelectric transducer (PCB 113B28) installed 50 mm above the nozzle exit. At each equivalence ratio condition, the combustion experiment is conducted for 8 s. The pressure signal is recorded with sampling rate 25 000 Hz, which is much faster than the mode of interest described below.

Figure 3.

(a) Schematic diagram of the annular combustor identical to that in Guk, Seo & Lee (Reference Guk, Seo and Lee2023). (b) Mode of interest featuring transverse acoustic oscillation, and (c) mean oscillation amplitude of this mode at varying methane–air equivalence ratio ($\phi$). (d,e) Bandpass-filtered pressure signal (black line) and its amplitude (blue line), and (f,g) p.d.f.s of the amplitude at (d,f) the fixed-point regime ($\phi =0.74$) and (e,g) the limit-cycle regime ($\phi =0.88$). MFC: mass flow controller.

By inspecting the oscillatory dynamics of the pressure signal, we found a distinct transverse mode at 355 Hz, which matches the result of the acoustic numerical simulation (figure 3(b), simulated using COMSOL Multiphysics v6.1). At this mode, a gradual increase in pressure oscillation amplitude is observed (figure 3c), implying a weakly nonlinear oscillation near a supercritical Hopf bifurcation. For diagnosing the thermoacoustic system both before and after the Hopf bifurcation, we select two experimental conditions, one in the fixed-point regime ($\phi =0.74$) and the other in the limit-cycle regime ($\phi =0.88$). The pressure time series and p.d.f.s of the oscillation amplitude at these conditions are shown in figures 3(d–g). Although the system identification using the stochastic Van der Pol equation (2.1) and the corresponding Fokker–Planck equation (2.4) often requires external stochastic forcing for inducing the NID (Lee et al. Reference Lee, Zhu, Li and Gupta2019, Reference Lee, Guan, Gupta and Li2020), we do not excite the system with additional noise, recognizing the strong inherent turbulence within the combustor (Noiray & Schuermans Reference Noiray and Schuermans2013; Lee et al. Reference Lee, Kim, Gupta and Li2021).

We employed time segmentation for the experimental data, which significantly improves training stability by increasing the number of samples at each time point. To be more precise, we partitioned each experimental dataset into eight time segments, assuming that the data are collected at the initial time point of each segment, i.e. $t=0, 1, \ldots, 7$. This assumption is substantiated by the fact that all experimental data originated from the stationary regime. It is worth mentioning that we have used an identical network for both the synthetic and experimental validations, differing only in the dataset used to train the log-likelihood loss function. In other words, we employed identical initial networks and maintained consistent loss functions, adjusting only the data for each specific problem.