A PINN approach for identifying governing parameters of noisy thermoacoustic systems

Abstract

Identifying the governing parameters of self-sustained oscillation is crucial for the diagnosis, prediction and control of thermoacoustic instabilities. In this paper, we propose and validate a novel method for computing the parameters of thermoacoustic oscillation in a stochastic environment, which exploits a physics-informed neural network (PINN). Specifically, we introduce a negative log-likelihood loss function that integrates the stochastic samples and the solution of the Fokker–Planck equation. The proposed framework is validated using the numerically generated signal and the experimental data obtained from an annular combustor, both before and after the supercritical Hopf bifurcation. The results of PINN-based system identification show good agreement with the actual system parameters and the original stochastic signal, with improved accuracy compared to established methods. To the best of our knowledge, this study constitutes the first demonstration of the PINN-inverse approach that uses the noise-induced dynamics of thermoacoustic systems, opening up new pathways for diagnosing and predicting the thermoacoustic behaviour of various combustion systems.

1. Introduction

1.1. Thermoacoustic oscillations via Hopf bifurcation

Despite extensive research spanning many decades, thermoacoustic oscillations remain a critical issue in developing and operating combustion systems such as rocket and gas turbine combustors. These oscillations occur from the positive feedback loop between the heat-release-rate fluctuations of an unsteady flame and the acoustic oscillations within the combustor (Lieuwen & Yang 2005). When these two types of oscillations align in phase, self-sustained thermoacoustic instability arises at one or more of the natural acoustic frequencies of the combustor via the Rayleigh mechanism (Magri, Juniper & Moeck 2020). In most cases, thermoacoustic oscillation is considered undesirable, as it exerts thermal stress on the combustion devices, potentially leading to mechanical damage or failure. Therefore, diagnosing the combustion system in terms of thermoacoustic instability is crucial in combustion systems so as to control the oscillation and prevent the pressure surge (Juniper & Sujith 2018; Krishnan et al. 2021).

Phenomenologically, birth and extinction of a thermoacoustic oscillation is often characterized as a Hopf bifurcation occurring between the fixed point and the limit cycle (Noiray & Schuermans 2013; Subramanian, Sujith & Wahi 2013). The normal form of the Hopf bifurcation in a nonlinear oscillating system is

(1.1)\begin{equation} \frac{\mathrm{d}{a}}{\mathrm{d}t}=k_1{a}+k_2a\,{|a|}^2, \end{equation}

where $a=|a|\,{\rm e}^{{\rm i}\phi }$ is the amplitude of oscillation, $t$ is time, $\phi$ is phase, and $k_1$ and $k_2$ are linear growth rate and nonlinear parameter, respectively. Negative and positive values of $k_1$ indicate fixed-point and limit-cycle regimes, respectively, with the Hopf bifurcation occurring at $k_1=0$ (Hopf point). Also, the negative nonlinear parameter ($k_2<0$) results in a supercritical Hopf bifurcation, where the limit cycle is observed only after the Hopf point ($k_1>0$). On the contrary, when $k_2$ is positive, subcritical Hopf bifurcation occurs, and the limit cycle can take place even at negative $k_1$. One or more higher-order nonlinear terms may be required to stabilize the oscillation in the latter situation. It is worth mentioning that in the hydrodynamics community, (1.1) is known as the Stuart–Landau equation (Landau 1944; Stuart 1960) characterizing the hydrodynamic instabilities in various types of flows, e.g. low-density jets (Raghu & Monkewitz 1991).

1.2. System identification using noise-induced dynamics for the diagnosis of thermoacoustic oscillations

When the amplitude oscillation in a thermoacoustic system is modelled with (1.1), it is crucial to identify the linear and nonlinear system parameters. In particular, the system parameters succinctly express the thermoacoustic dynamics, thus the system can be diagnosed using these values. The computation of system parameters can be conducted via system identification, which involves statistical analysis of the system's output signal. An established method of system identification near a Hopf bifurcation uses the noise-induced dynamics (NID) of the system (Lee et al. 2019; Kabiraj, Vishnoi & Saurabh 2020). Opposing the conventional view of regarding noise as a signal contamination, researchers have shown that a stochastically perturbed system can reveal more information about its original deterministic dynamics than a noise-free system via the NID (Pikovsky & Kurths 1997; Ushakov et al. 2005; Kabiraj et al. 2015). Therefore, in the system identification framework mentioned above, a stochastic differential equation (SDE) is postulated to capture the noise-induced behaviour of the system. Near the Hopf point where the growth rate is small, the system is weakly nonlinear, and the stochastic averaging can be applied to this SDE, yielding an equivalent Fokker–Planck equation describing the drift and diffusion dynamics of the thermoacoustic system. Siegert, Friedrich & Peinke (1998) have proposed that drift and diffusion terms of the Fokker–Planck equation can be extracted by analysing the time correlation of the noisy data, establishing an output-only system identification method. Noiray & Schuermans (2013) first applied this framework to compute linear and nonlinear parameters of the stochastic Van der Pol equation governing thermoacoustic oscillations in a gas-turbine combustor. Noiray & Denisov (2017) later validated this stochastic system identification framework by applying periodic feedback control to the combustor at both the stable and unstable regimes. Consequently, the system identification method using the stochastic Van der Pol equation and the corresponding Fokker–Planck equation is applied to numerous combustion systems exhibiting thermoacoustic oscillation, including both laminar and turbulent combustors (Boujo & Noiray 2017; Bonciolini, Boujo & Noiray 2017; Lee et al. 2020, 2021).

Despite the robustness of the NID-based system identification framework proven in various thermoacoustic systems, there exists an intrinsic limitation that hinders the accuracy of the extracted parameters. First, established methods of system identification that exploit the NID of the system are based on the probability mass function, counting discrete samples of the time-shifted signal. This causes the inherent discrete approximation error, especially when the input sample size is small. Furthermore, while computing the time correlation of the signal, the noise is assumed to be perfectly memoryless. In practical thermoacoustic systems, however, Markovian assumption can be severely impaired due to the coarse sampling rate and the signal filtering. In such situations, adverse finite-time effects arise, and further adjustment of the computed parameters may be required, specifically by incorporating the optimization scheme based on adjoint equations (Lade 2009; Boujo & Noiray 2017; Lee, Kim & Park 2023b). To the best of our knowledge, a NID-based system identification method that is unaffected by the above-mentioned limitations due to discrete time-shift computation has yet to be developed.

1.3. Physics-informed neural networks for the system identification of thermoacoustic systems

Along with recent advances in machine learning, neural networks are increasingly applied for the analysis of thermoacoustic systems. Selimefendigil & Polifke (2011) developed a low-order model in the frequency domain for predicting thermoacoustic limit cycles using the feed-forward neural network identification method. Consequently, neural networks are exploited for deducing heat release models (Jaensch & Polifke 2017) and nonlinear flame responses (Tathawadekar et al. 2021) in thermoacoustic systems. Recently, Nóvoa & Magri (2022) used an echo state network, a reservoir-computing-based recurrent neural network, for the real-time bias-aware estimation of the states and parameters of a numerical Rijke tube model. Nóvoa, Racca & Magri (2024) later generalized this echo state network with a regularized bias-aware ensemble Kalman filter. Building on these established methods, we aim to develop a physics-informed neural network (PINN) for solving the NID-based inverse problem introduced in § 1.2.

A PINN is a neural network that solves forward and inverse problems while respecting the governing equations in the form of partial differential equations (PDEs) (Lagaris, Likas & Fotiadis 1998; Karniadakis et al. 2021). By reason of their ability to provide desirable solutions to ill-posed problems, PINNs have been used widely for the inverse modelling of dynamical systems. For example, Chen et al. (2020) employed a PINN to tackle ill-posed inverse scattering problems in nano-optics photonic metamaterials. Jagtap et al. (2022) used the extended PINN framework for addressing inverse supersonic compressible flow problems. For enhancing the training accuracy through inductive bias, hard-constrained PINNs for inverse problems are proposed in Lu et al. (2021). Recently, Ozan & Magri (2023) integrated a hard-constrained PINN and the Galerkin decomposition technique to model nonlinear acoustics in a prototypical thermoacoustic system. One can find further theoretical foundations for the PINN approach to inverse problems in Mishra & Molinaro (2022) and Zhang, Li & Liu (2023).

In order to model the noise-induced behaviour of the thermoacoustic system, it is essential to incorporate SDEs into the PINN framework. Although not applied to thermoacoustic systems, recent studies have suggested that PINNs can be adopted to solve inverse problems of the SDEs (Xu & Darve 2021; Shin & Choi 2023), opening up possibilities for system identification in noise-perturbed combustors. However, only a limited number of studies have integrated stochastic samples with averaged PDEs, which is essential for system identification using the NID. For instance, Chen et al. (2021) proposed a PINN approach for solving inverse problems involving NID through the Fokker–Planck equation using discrete particle observations. The authors employed the Kullback–Leibler divergence from the observed empirical distribution to the neural network solution as a loss function. We, on the other hand, seek to minimize the discrepancy between the discrete samples and the analytical Fokker–Planck equation using the maximum likelihood approach, aiming to obtain more stable solutions suitable for system identification in noisy combustors.

1.4. Contributions of the present study

In this paper, we aim to develop and validate an NID-based system identification method that does not require the computation of the time correlation of the signal, and thus does not suffer from the discrete approximation error and the finite-time effect. Specifically, we aim to exploit a PINN linking the combustor signal and the Fokker–Planck equation equivalent to the stochastic Van der Pol equation for solving the inverse problem (i.e. system identification). A key mathematical distinction in our framework is the use of the maximum likelihood approach to integrate discrete samples with the Fokker–Planck equation, directly incorporating stochastic samples from the original SDE.

Below, we present our mathematical framework in § 2, focusing on the system model and PINN design. Next, we describe the numerical and experimental data for validating the proposed method in § 3. We then show the system identification results and the following discussion in § 4, before concluding in § 5.

2. Mathematical framework

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. (2023b) 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 1981; Zhu & Yu 1987). 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. 2021; Indlekofer et al. 2022). 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. 2021; Xu & Darve 2021; Shin & Choi 2023) 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 2023). 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 2014) 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. Data for PINN validation

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 (2023).

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. (2023) (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 (2023). (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. 2019, 2020), we do not excite the system with additional noise, recognizing the strong inherent turbulence within the combustor (Noiray & Schuermans 2013; Lee et al. 2021).

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.

4. Results and discussion

In this section, we assess the results of PINN-based system identification in terms of the parameter accuracy and the likelihood of the reconstructed p.d.f. We also present comparisons of the results to the existing methods of system identification that use the NID, which use extrapolation of the time correlation (Lade 2009; Noiray & Schuermans 2013), and the same method with adjoint-based optimization (Boujo & Noiray 2017). These methods will be denoted SI-ext and SI-opt, respectively, in the following figures. In both of the established methods, drift and diffusion terms are estimated using the equation

(4.1)\begin{equation} D^{(n)}(a) = \lim_{\tau \rightarrow 0} D_{\tau}^{(n)}(a), \end{equation}

where $D_{\tau }^{(n)}(a)$ is defined by

(4.2)\begin{equation} D_{\tau}^{(n)}(a) = \frac{1}{n!\,\tau} \int_0^{\infty} (A-a)^n\,P(A, t + \tau \mid a, t) \,{\rm d}A, \end{equation}

where $P(A, t+\tau \mid a,t)$ is the conditional p.d.f. of the pressure oscillation amplitude being $A$ at time $t + \tau$ given $a$ at time $t$. In the adjoint-based optimization framework (SI-opt), the discrepancy between $D_{\tau }^{(n)}(a)$ computed from the obtained parameters and the experimental measurement is minimized using the optimization scheme (Boujo & Noiray 2017). It should be noted that neither the time-shifted drift/diffusion term estimation nor the adjoint-based optimization is applied in the proposed PINN framework, and these methods are used only for comparing the extracted system parameters in this paper. For further information about these existing methods, readers may refer to Lee et al. (2023b).

First, for identifying the governing parameters of the synthetic limit-cycle data, we initialized $\epsilon$, $\alpha$ and $d$ to be zero and applied the Xavier initialization for the network parameters before training, following the guidelines outlined by Glorot & Bengio (2010). For the computation of loss functions $\mathcal {L}_R$ and $\mathcal {L}_{mass}$, we sampled 10 000 collocation points from the domain $[0, A] \times [0, T]$. Additionally, we sampled 100 points for each component of $\mathcal {L}_{BC}$. We set the initial learning rate at $10^{-4}$, and reduced it systematically using a learning rate scheduler throughout the training process.

Figures 4(a–f) show the evolution of the parameters as training iterates when using synthetic transient data. While the parameters $\epsilon,\alpha,d$ converge rapidly in the limit-cycle cases ($\epsilon >0$), fixed-point cases ($\epsilon <0$) tend to require greater number of iterations for convergence, especially when the system deviates further from the Hopf point ($\epsilon =0$). This is because the fixed-point data are generally less deterministic compared to limit-cycle counterparts, and thus have a wider choice of system parameters for minimizing the loss function. Such a trend is also found in existing system identification methods, as shown in table 1. Nevertheless, in all cases, PINN-based system identification shows substantial improvement from the existing system identification methods in terms of the average relative error of the identified parameters (see table 1). As a result, we were able to reconstruct accurately the time evolution of the oscillation amplitude and its corresponding p.d.f. for both the fixed-point and limit-cycle data, as shown in figures 4(g–r). Considering that just ten stochastic signals (figure 2) are used as the input data in each case, the accuracy of the present system identification framework is recognizable. Again, the parameters identified from the PINN approach are closer to the true values in all test cases, as depicted in figure 5.

Figure 4. (a–f) The convergence history of $\epsilon, \alpha, d$ during PINN-based system identification. (g–l) Time evolution of the oscillation amplitude obtained from 500 numerical repetitions using the identified parameters, and (m–r) the corresponding p.d.f.s. The results are obtained from six sets of synthetic data: (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. Grey dashed lines indicate true values of $\epsilon$, $\alpha$ and $d$.

Figure 5. Parameters of the stochastic Van der Pol equation computed from extrapolation-based system identification (SI-ext, blue line with circular markers), system identification with adjoint-based optimization (SI-opt, green line with cross markers), and the present method (PINN, red line with diamond markers).

Table 1. Parameters identified from synthetic data using the following methods: extrapolation-based system identification (SI-ext), adjoint-optimization-based system identification (SI-opt) and PINN. Here, $E$ is the average relative error of the identified parameters ($\epsilon,\alpha,d$) compared to true values.

Next, we present the PINN-based system identification results for the annular combustor data in table 2 and figure 6. Because the true values of the governing parameters of the thermoacoustic oscillation are unknown, we alternatively assess the average likelihood values $L_{avg}=-({1}/{M})\sum _{i,j=1}^{N,M}\log P(a_i, t_j)$ of the analytical solution computed from the identified parameters. It can be seen from table 2 that $L_{avg}$ is greatest in the PINN-based system identification for both the fixed-point and limit-cycle data. This indicates that the parameters computed from the PINN approach best represent the experimental p.d.f., especially in the limit-cycle regime (see figure 6 for comparison).

Figure 6. Comparison between the experimental data (grey bars) and stationary analytical Fokker–Planck solutions using parameters computed from extrapolation-based system identification (SI-ext, blue dashed line), system identification with adjoint-based optimization (SI-opt, green dotted line), and the present method (PINN, red continuous line) at (a) the fixed-point regime ($\phi =0.74$) and (b) the limit-cycle regime ($\phi =0.88$). The stationary analytical solution could not be computed with SI-ext parameters at $\phi =0.74$ because of the diverging solution ($\alpha >0$). Here, $\phi$ is the methane–air equivalence ratio.

Table 2. Parameters identified from experimental data at the equivalence ratios (a) $\phi =0.74$ and (b) $\phi =0.88$, using the following methods: extrapolation-based system identification (SI-ext), adjoint-optimization-based system identification (SI-opt) and PINN. Here, $L_{avg}$ is the likelihood of the p.d.f. reconstructed from identified parameters $\epsilon,\alpha,d$.

Finally, we reconstruct the phase portrait of the limit cycle from the computed parameters of thermoacoustic oscillation. We use the delay embedding technique proposed by Takens (1981), which enables the reconstruction of the phase portrait using just a single time series with a time delay ($\tau$). This method is used widely to examine the dynamics of thermoacoustic oscillations (Kashinath, Li & Juniper 2018; Lee et al. 2020). Following Fraser & Swinney (1986), we choose the minimum $\tau$ that minimizes the average mutual information. The reconstructed phase portraits shown in figure 7 reveal that the parameters identified from the PINN-based framework best capture the dynamics of the experimental data in both the fixed-point and limit-cycle regimes, reconfirming the robustness of the present method.

Figure 7. Phase portraits obtained from (a–c,e–g) numerical simulations and (d,h) experimental data. Numerical results are computed using the system identification results of the experimental data obtained from (a,e) extrapolation-based SI (SI-ext), (b,f) SI with adjoint-based optimization (SI-opt), and (c,g) the present method (PINN). Grey scatter dots are obtained from the stochastic signal, while blue lines are computed from deterministic ($d=0$) simulations. Cross markers in (b,c) indicate analytical fixed points, while the fixed point could not be obtained with SI-ext parameters at $\phi =0.74$ because of the diverging solution ($\alpha >0$). Here, $\tau$ is the time delay computed from the minimum average mutual information (Fraser & Swinney 1986), and $\phi$ is the methane–air equivalence ratio.

5. Conclusions

In this study, we performed NID-based system identification of a thermoacoustic oscillator in fixed-point and limit-cycle regimes using the PINN approach. From numerical and experimental validation, we found that the system identification using PINN leads to better computation performance than existing system identification methods that incorporate NID, in terms of parameter accuracy and likelihood. This is the first time that PINN has been applied for diagnosing thermoacoustic oscillations modelled with a stochastic oscillator equation, to the best of our knowledge. A major implication of this study is that the proposed framework could identify accurately the parameters of the stochastic self-sustained oscillation using the output-only method without requiring information about the input signal or the adjoint-based optimization schemes. Thus an efficient and versatile system identification could be performed in a noisy thermoacoustic system by capturing the NID from the PINN framework.

There are three categories of studies that can be explored in the future. First, albeit not demonstrated in this paper, parameters of the stochastic Van der Pol equation identified from the present framework can be used to predict and control the thermoacoustic instabilities. Specifically, the Hopf point and the post-bifurcation dynamics can be predicted via the extrapolation of the identified parameters (Lee et al. 2020), and the feedback control can be conducted using these parameters. Second, the present study dealt with the occurrence of thermoacoustic oscillation via the supercritical Hopf bifurcation. In future research, other routes to instabilities involving higher-order dynamics, such as subcritical Hopf bifurcation (Gopalakrishnan et al. 2016), intermittency (Nair, Thampi & Sujith 2014) and period doubling (Subramanian et al. 2010), can be studied using the present PINN-based method. Finally, considering that the system model used in this study is purely phenomenological, the present framework can be applied to other physical systems exhibiting stochastic self-sustained oscillations. For example, one can apply PINN-based system identification to diagnose hydrodynamically oscillating low-density jets (Zhu, Gupta & Li 2017; Lee et al. 2019; Park & Lee 2024) or electrically oscillating plasma in Hall-effect thrusters (Han et al. 2023; Lee et al. 2023a).