A real-time digital twin of azimuthal thermoacoustic instabilities

Abstract

When they occur, azimuthal thermoacoustic oscillations can detrimentally affect the safe operation of gas turbines and aeroengines. We develop a real-time digital twin of azimuthal thermoacoustics of a hydrogen-based annular combustor. The digital twin seamlessly combines two sources of information about the system: (i) a physics-based low-order model; and (ii) raw and sparse experimental data from microphones, which contain both aleatoric noise and turbulent fluctuations. First, we derive a low-order thermoacoustic model for azimuthal instabilities, which is deterministic. Second, we propose a real-time data assimilation framework to infer the acoustic pressure, the physical parameters, and the model bias and measurement shift simultaneously. This is the bias-regularized ensemble Kalman filter, for which we find an analytical solution that solves the optimization problem. Third, we propose a reservoir computer, which infers both the model bias and measurement shift to close the assimilation equations. Fourth, we propose a real-time digital twin of the azimuthal thermoacoustic dynamics of a laboratory hydrogen-based annular combustor for a variety of equivalence ratios. We find that the real-time digital twin (i) autonomously predicts azimuthal dynamics, in contrast to bias-unregularized methods; (ii) uncovers the physical acoustic pressure from the raw data, i.e. it acts as a physics-based filter; (iii) is a time-varying parameter system, which generalizes existing models that have constant parameters, and capture only slow-varying variables. The digital twin generalizes to all equivalence ratios, which bridges the gap of existing models. This work opens new opportunities for real-time digital twinning of multi-physics problems.

1. Introduction

Thermoacoustic instabilities are a multi-physics phenomenon, which is caused by the constructive coupling between hydrodynamics, unsteady heat released by flames and acoustics (e.g. Culick 1988; Paschereit, Gutmark & Weisenstein 1998; Lieuwen et al. 2001; Candel et al. 2009; Poinsot 2017; Juniper & Sujith 2018; Magri, Schmid & Moeck 2023; Silva 2023). A thermoacoustic instability can arise when the heat released by the flames is sufficiently in phase with the acoustic pressure (Rayleigh 1878; Magri, Juniper & Moeck 2020). If uncontrolled or not prevented, instabilities grow into large-amplitude pressure oscillations, which can detrimentally affect gas turbine operating regimes, cause structural damage and fatigue and, in the worst-case scenario, shake the engine and its components apart (e.g. Candel 2002; Dowling & Morgans 2005; Culick 2006; Lieuwen 2012). In aeroengines the flame holders are arranged in annular configurations to increase the power density (e.g. Krebs et al. 2002). Despite these annular configurations being nominally rotationally symmetric, large-amplitude azimuthal thermoacoustic oscillations can spontaneously occur and break dynamical symmetry (e.g. Morgans & Stow 2007; Noiray, Bothien & Schuermans 2011; Noiray & Schuermans 2013b; Bauerheim, Nicoud & Poinsot 2016; Humbert et al. 2023a; Magri et al. 2023). Thermoacoustic instabilities in annular combustors have intricate dynamics, which can be grouped into (e.g. Magri et al. 2023) (a) spinning, if the nodal lines rotate azimuthally (typical of rotationally symmetric configurations); (b) standing, if the nodal lines have, on average, a fixed orientation (typical of rotationally asymmetric configurations); and (c) mixed, if the nodal lines switch between the two former states (typical of weakly asymmetric configurations) (e.g. Schuermans, Paschereit & Monkewitz 2006; Noiray et al. 2011; Worth & Dawson 2013a). Azimuthal instabilities have been investigated by high-fidelity simulations (e.g. Wolf et al. 2012), by experimental campaigns in atmospheric and pressurized rigs (e.g. Bourgouin et al. 2013; Worth & Dawson 2013b; Mazur et al. 2019, 2021; Ahn et al. 2022; Indlekofer et al. 2022) and heavy-duty gas turbines (Noiray & Schuermans 2013b), and by theoretical studies (e.g. Moeck, Paul & Paschereit 2010; Noiray et al. 2011; Ghirardo & Juniper 2013; Duran & Morgans 2015; Bauerheim et al. 2016; Laera et al. 2017; Mensah et al. 2019; Murthy et al. 2019; Faure-Beaulieu & Noiray 2020). Because the dynamics of azimuthal thermoacoustics are not yet fully understood (Aguilar Pérez et al. 2021; Faure-Beaulieu et al. 2021b), the understanding, modelling and control of azimuthal oscillations is an active area of research.

Experimental campaigns were performed to gain insight into the physical mechanisms and behaviour of azimuthal instabilities in a prototypical annular combustor with electrically heated gauzes (Moeck et al. 2010). In a model annular gas turbine combustor with ethylene–air flames, Worth & Dawson (2013b) investigated the flame dynamics and heat release, and how it coupled with the acoustics (O'Connor, Worth & Dawson 2013; Worth & Dawson 2013a). They found that varying the burner spacing, to promote or suppress flame–flame interactions, resulted in changes to the amplitude and frequency of the azimuthal modes. They also varied the burner swirl directions, which impacted the preferred mode selection. The spontaneous symmetry breaking of thermoacoustic eigenmodes was experimentally analysed in Indlekofer et al. (2022), who found that small imperfections in the rotational symmetry of the annular combustor were magnified at the supercritical Hopf bifurcation point that separated resonant from limit cycle oscillations in the form of a standing mode oriented at an azimuthal angle defined by the rotational asymmetry. The nonlinear dynamics of beating modes was experimentally discovered and modelled by Faure-Beaulieu et al. (2021b). The authors found that for some combinations of small asymmetries of the resistive and reactive components of the thermoacoustic system, purely spinning limit cycles became unstable, and heteroclinic orbits between the corresponding saddle points led to beating oscillations with periodic changes in the spin direction. The previous works focused on statistically stationary regimes. The analysis of slowly varying operating conditions were investigated in Indlekofer et al. (2021), who unravelled dynamic hystereses of the thermoacoustic state of the system. Other experimental and theoretical investigations focused on analysing the intermittent behaviour (Faure-Beaulieu et al. 2021a; Roy et al. 2021). In the latter work, the solutions of the Fokker–Plank equation, which governs the probability of observing an instantaneous limit cycle state dominated by a standing or spinning mode, enabled the prediction of the first passage time statistics between erratic changes in spin direction.

The intricate linear and nonlinear dynamics of azimuthal thermoacoustic oscillations spurred interest in physics-based low-order modelling and control (e.g. Morgans & Stow 2007; Illingworth & Morgans 2010; Humbert et al. 2023b). The azimuthal dynamics can be qualitatively described by a one-dimensional wave-like equation (Noiray et al. 2011; Ghirardo & Juniper 2013; Bauerheim et al. 2014; Bothien, Noiray & Schuermans 2015; Yang, Laera & Morgans 2019), which can also include a stochastic forcing term to model the effect of the turbulent fluctuations and noise (Noiray & Schuermans 2013a; Orchini et al. 2020). In the frequency domain, azimuthal oscillations were investigated with eigenvalue sensitivity (Magri, Bauerheim & Juniper 2016a; Magri et al. 2016b; Mensah et al. 2019; Orchini et al. 2020). These works concluded that traditional eigenvalue sensitivity analysis needed to be extended to tackle degenerate pairs of azimuthal modes, as reviewed in Magri et al. (2023). A review of azimuthal thermoacoustic modelling can be found in Bauerheim et al. (2016). In the time domain, the typical approach is to develop models that describe the acoustic state, which is also referred to as the slow-varying variable approach based on quaternions (e.g. Ghirardo & Bothien 2018) or as the generalized Bloch sphere representation (Magri et al. 2023). In this approach, all fast-varying dynamics, such as aleatoric noise and turbulent fluctuations, are filtered out of the data and modelled as a stochastic forcing term. The equations model the amplitude of the acoustic pressure envelope, the temporal phase drift and the angles defining rotational and reflectional symmetry breaking (Ghirardo & Bothien 2018). Faure-Beaulieu & Noiray (2020) introduced the slow-varying variables into an stochastic wave equation, which was averaged in space and time to obtain coupled Langevin equations. The low-order model parameters were estimated via Langevin regression, which is a regression method used in turbulent environments (e.g. Siegert, Friedrich & Peinke 1998; Noiray 2017; Boujo et al. 2020; Callaham et al. 2021). Slow-varying variable models provide qualitatively accurate representations of azimuthal oscillations. However, they are not suitable for real-time applications, in which we need models that can receive data from sensors as raw inputs. Therefore, this work focuses on models describing the fast-varying quantities, which model the evolution of the modal amplitudes through coupled Van der Pol oscillators (e.g. Noiray et al. 2011). Nonetheless, the discussed low-order modelling methods are offline, i.e. they infer the model parameters from the data in a post-processing stage and they identify one parameter at a time (i.e. they do not estimate all the parameters in one computation). This means that the physical parameters of the literature are not necessarily optimal.

On the one hand, physics-based low-order models are qualitatively accurate, but they are quantitatively inaccurate due to modelling assumptions and approximations (Magri & Doan 2020). The low-order models’ state and parameters are affected by both aleatoric uncertainties and epistemic errors, i.e. model biases (Nóvoa, Racca & Magri 2024). On the other hand, experimental data can provide reliable information about the system, but they are typically noisy and sparse (e.g. Magri & Doan 2020). Experimental measurements can be affected by aleatoric uncertainties due to environmental and instrumental noise, and systematic errors in the sensors. To rigorously combine the two sources of information (low-order modelling and experimental data) to improve the knowledge on the system, data assimilation (DA) comes into play (e.g. Tarantola 2005; Evensen 2009). Data assimilation for thermoacoustics was introduced with a variational approach (offline) by Traverso & Magri (2019) with a real-time (or sequential) approach by Nóvoa & Magri (2022), who deployed an ensemble square-root Kalman filter to infer the pressure and physical parameters. Sequential DA has also been successfully applied to reacting flows (e.g. Labahn et al. 2019; Yu et al. 2019a; Donato, Galletti & Parente 2024), turbulence modelling (e.g. Colburn, Cessna & Bewley 2011; Majda & Harlim 2012; Gao et al. 2017; Magri & Doan 2020; Hansen, Brouzet & Ihme 2024) and acoustics (e.g. Kolouri, Azimi-Sadjadi & Ziemann 2013; Wang et al. 2021), among others. Apart from Nóvoa & Magri (2022), these works implemented classical DA formulations, which assume that the uncertainties are aleatoric and unbiased (e.g. Dee 2005; Laloyaux et al. 2020). However, as shown in Nóvoa et al. (2024), a DA framework provides an optimal state and set of parameters when the model biases are modelled. Model bias estimation in real-time DA is traditionally based on the separate-bias Kalman filter scheme (Friedland 1969). This bias-aware filter augments the dynamical system with a parameterized model for the bias, and then solves two state and parameter estimation problems: one for the physical model and another for the bias model (e.g. Ignagni 1990; Dee & da Silva 1998; da Silva & Colonius 2020). However, the separate-bias Kalman filter relies on the knowledge of the bias functional form a priori. Also, as explained in Nóvoa et al. (2024), the filter is unregularized, which can lead to unrealistically large estimates of the model bias and does not ensure that the bias is unique. In other words, the separate-bias Kalman filter may be ill-posed. To overcome these limitations, Nóvoa et al. (2024) proposed a regularized bias-aware ensemble data assimilation framework, for which they derived an analytical solution: the bias-regularized ensemble Kalman filter (r-EnKF). The r-EnKF is a real-time DA method that not only accounts for the model bias, but also regularizes the norm of the bias. The regularization makes the algorithm stable and the solution of the problem unique. In the limit of a perfect model (i.e. zero model bias), the r-EnKF recovers the classical bias-unaware ensemble Kalman filter (EnKF) introduced by Evensen (2003). The model bias was inferred by an echo state network (ESN), which is a universal approximator of time-varying functions (e.g. Grigoryeva & Ortega 2018) that can adaptively estimate the model bias with no major assumption regarding the functional form (Nóvoa et al. 2024). Echo state networks are suitable for real-time digital twins because their training consists of solving a linear system, which is computationally inexpensive in contrast to the back-propagation methods required in other machine learning methods (e.g. Bonavita & Laloyaux 2020; Brajard et al. 2020). The r-EnKF of Nóvoa et al. (2024) assumed that the measurements were unbiased, which is an assumption that we relax in this paper to assimilate real experimental data.

1.1. Objectives and structure

The objective of this work is to develop a real-time digital twin of azimuthal thermoacoustic instabilities. Figure 3 pictorially describes the proposed digital twin framework. Real-time digital twins are adaptive models, which are designed to predict the behaviour of their physical counterpart by assimilating data from sensors when they become available. We need four ingredients to design a real-time digital twin (i) data from sensors (from noisy and sparse microphones), (ii) a qualitative low-order model, which should capture the fast-varying acoustic dynamics rather than the slow-varying states, (iii) an estimator of both the model bias and measurement shift, (iv) a statistical DA method that optimally combines the data and the model on the fly (in real time) to infer the physical states and parameters. Section 2 describes the available experimental data, as well as the physics-based low-order model used to describe the thermoacoustic problem. The real-time assimilation framework is detailed in § 3. The proposed ESN is described in § 4. Section 6 shows the results of the real-time digital twin, and compares the performance of the bias-unregularized EnKF with the r-EnKF. § 7 ends the paper.

2. Azimuthal thermoacoustic instabilities

We model in real time azimuthal thermoacoustic instabilities in hydrogen-based annular combustors. In this section we describe both the experimental data (§ 2.1) and the physics-based low-order model (§ 2.2).

2.1. Experimental set-up and data

We use the experimental data of Faure-Beaulieu et al. (2021a) and Indlekofer et al. (2022). The experimental set-up (figure 1) consists of an annular combustor with 12 equally spaced burners with premixed flames, which are fuelled with a mixture of 70/30 % ${\rm H}_2/{\rm CH}_4$ by power. The operating conditions are atmospheric and the thermal power is fixed at 72 kW. The equivalence ratios considered are $\varPhi =\{0.4875, 0.5125, 0.5375, 0.5625\}$. When $\varPhi >0.5$, the system is thermoacoustically unstable with self-sustained oscillations that peak at approximately 1090 Hz (Indlekofer et al. 2022). The acoustic pressure data was recorded at a sampling rate of 51.2 kHz by Kulite pressure transducers (XCS-093-05D) at four azimuthal locations $\theta =\{0^\circ, 60^\circ, 120^\circ, 240^\circ \}$.

Figure 1. Side and top views of the experimental set-up of the annular combustor. Adapted from Faure-Beaulieu et al. (2021b).

Figure 2 shows the time series, histograms and power spectral density (PSD) of the raw and post-processed experimental data at $\varPhi =0.5125$. The PSD of the post-processed data is obtained after applying a Butterworth band-pass filter to the raw data around the frequency of instability, which approximately isolates the frequencies between 1050 and 1150 Hz (figure 2 III). The pressure measurements have a mean value that is different from zero because of large-scale flow structures, which are not correlated to the thermoacoustic dynamics. The non-zero mean pressure is referred to as ‘measurement shift’, as further explained in § 3.1. In a digital twin that is coupled with sensors’ measurements in real time, the data are the raw pressure signals that are assimilated into the low-order model (figure 2a). Because the raw pressure contains aleatoric noise, measurement shift and turbulent flow fluctuations, we refer to the post-processed data (figure 2b) as the ‘presumed ground truth’ or ‘presumed acoustic state’, that is, the physical acoustic state that is the presumed target prediction.

Figure 2. Pressure at the four measurement locations with equivalence ratio $\varPhi =0.5125$ (thermoacoustically unstable configuration). (a) Raw experimental data, i.e. the observables, and (b) post-processed data, which are briefly referred to as the presumed truth. (I) Time series of the fast-varying oscillations, (II) histogram of the acoustic pressure in the time window for DA, $t\in [0.5, 0.85]$, and (III) comparison of the power spectral density of the raw data (light colour) and the presumed truth (dark colour). The horizontal lines in (II) indicate the mean of the histograms. The thermoacoustic limit cycle is a standing acoustic mode with a nodal line near $\theta =120^\circ$ (mono-modal histogram).

Figure 3. Schematic of the proposed digital twin framework. (a) The physical and digital systems evolve in time (left to right) independently. The digital system is composed of a physical model with uncertain state and parameters, and an estimator of the model bias and innovations (from which the measurement shift is estimated). When sensor data (crosses) become available, they are combined with the estimates from the digital twin using the r-EnKF to update the digital system. (b) Diagram of the r-EnKF update performed sequentially every time that data become available.

2.2. A qualitative and physics-based low-order model

The dynamics of the azimuthal acoustics can be qualitatively modelled by a one-dimensional wave-like equation (e.g. Faure-Beaulieu & Noiray 2020; Indlekofer et al. 2022)

(2.1)\begin{gather} \dfrac{\partial^2p}{\partial t^2} + \zeta\dfrac{\partial p}{\partial t} - \left[1+\epsilon\cos{\left(2(\theta-\varTheta_\epsilon)\right)}\right]\dfrac{c^2}{r^2}\dfrac{\partial^2p}{\partial \theta^2} = (\gamma - 1) \dfrac{\partial \dot{q}}{\partial t}, \end{gather}

(2.2)\begin{gather} \quad \text{with}\ (\gamma - 1) \dfrac{\partial \dot{q}}{\partial t} = \beta\left[1+c_2\cos{(2(\theta-\varTheta_\beta))}\right]p - \kappa p^3, \end{gather}

where $p$ is the acoustic pressure; $\zeta$ is the acoustic damping; $c$ is the speed of sound; $r$ is the mean radius of the annulus; $\gamma$ is the heat capacity ratio; $\epsilon$ and $\varTheta _\epsilon$ are the amplitude and phase of the reactive symmetry, respectively (Indlekofer et al. 2022); $\dot {q}$ is the coherent component of the fluctuations of the heat release rate (Noiray 2017), which is divided into a nonlinear cubic term, which models the saturation of the flame response weighted by the parameter $\kappa$; and a linear response to the acoustic perturbations, which is weighted by the heat release strength $\beta$, the resistive asymmetry intensity, $c_2$, and direction of maximum root-mean-square (r.m.s.) acoustic pressure, $\varTheta _\beta$. The flame response model (2.2) is an accurate approximation in the vicinity of a Hopf bifurcation (e.g. Lieuwen 2003; Noiray 2017). To further reduce the complexity of (2.2), we project the acoustic pressure on the degenerate pair of eigenmodes of the homogeneous wave equation (Noiray et al. 2011; Noiray & Schuermans 2013b; Magri et al. 2023)

(2.3)\begin{equation} p(t, \theta) = \eta_a(t) \cos\theta + \eta_b(t)\sin\theta, \end{equation}

where $\eta _a$ and $\eta _b$ are the acoustic velocity amplitudes. Substituting (2.3) into (2.1) and averaging yield the governing equations, which consist of a set of nonlinearly coupled oscillators

(2.4a)\begin{align} \ddot{\eta}_a &={-}\omega^2\left[\eta_a \left(1 + \frac{\epsilon}{2}\cos{(2 \varTheta_\epsilon)}\right) + \eta_b \frac{\epsilon}{2} \sin{(2\varTheta_\epsilon)}\right] \nonumber\\ &\quad +\dot{\eta}_a\left[2\nu+\frac{c_2\beta}{2}\cos{(2\varTheta_\beta)}-\frac{3\kappa}{4}(3\eta_a^2+\eta_b^2)\right]+ \dot{\eta}_b\left[\frac{c_2\beta}{2}\sin{(2\varTheta_\beta)}- \frac{3}{2}\kappa\eta_b\eta_a\right], \end{align}

(2.4b)\begin{align} \ddot{\eta}_b &={-}\omega^2\left[\eta_b \left(1 - \frac{\epsilon}{2}\cos{(2 \varTheta_\epsilon)}\right) + \eta_a \frac{\epsilon}{2} \sin{(2\varTheta_\epsilon)}\right] \nonumber\\ &\quad +\dot{\eta}_b\left[2\nu-\frac{c_2\beta}{2}\cos{(2\varTheta_\beta)}-\frac{3\kappa}{4}(3\eta_b^2+\eta_a^2)\right]+ \dot{\eta}_a\left[\frac{c_2\beta}{2}\sin{(2\varTheta_\beta)}- \frac{3}{2}\kappa\eta_b\eta_a\right], \end{align}

where $\nu =(\beta -\zeta )/2$ is the linear growth rate of the pressure amplitude in absence of asymmetries, i.e. when $c_2=0$ (Indlekofer et al. 2022). Equation (2.4) can be written in a compact notation with a nonlinear state-space formulation

(2.5) \begin{equation} \left.\begin{array}{c@{}} {\rm d}{\boldsymbol{\mathit{\phi}}}= \mathcal{F}\left(\boldsymbol{\mathit{\phi}}, \boldsymbol{\alpha} \right) {\rm d}t, \\ {\boldsymbol{\mathit{q}}} = \mathcal{M}({\boldsymbol{\mathit{\theta}}}, {\boldsymbol{\mathit{\phi}}}), \end{array} \right\} \end{equation}

where ${\boldsymbol {\mathit {\phi }}}=[\eta _a; \dot {\eta }_a; \eta _b; \dot {\eta }_b]\in \mathbb {R}^{N_\phi }$ is the state vector; $\mathcal {F}:\mathbb {R}^{N_\phi }\rightarrow \mathbb {R}^{N_\phi }$ is the nonlinear operator that represents (2.4) and ${\boldsymbol {\mathit {\alpha }}}=[\nu ; c_2\beta ; \omega ; \kappa ; \epsilon ; \varTheta _\beta ; \varTheta _\epsilon ]\in \mathbb {R}^{N_\alpha }$ are the system's parameters (the operator $[~;~]$ indicates vertical concatenation and we use column vectors throughout the paper). The vector ${\boldsymbol {\mathit {q}}}$ describes the model observables, which are the acoustic pressures at the azimuthal locations ${\boldsymbol {\mathit {\theta }}}= \{0^\circ, 60^\circ, 120^\circ, 240^\circ \}$. The model observables are computed through the measurement operator $\mathcal {M}:\mathbb {R}^{N_\phi }\rightarrow \mathbb {R}^{N_q}$, which maps the state variables into the observable space through (2.3).

3. Real-time DA

Real-time DA makes qualitative models quantitatively accurate every time that sensors’ measurements (data) become available (Magri & Doan 2020). At a measurement time $t_k$, the model's state, ${\boldsymbol {\mathit {\phi }}}_k$, and parameters, ${\boldsymbol {\mathit {\alpha }}}_k$, are inferred by combining two sources of information, that is, the measurement, ${\boldsymbol {\mathit {d}}}_k$, and the model forecast, ${\boldsymbol {\mathit {q}}}_k$. A robust inference process must also filter out both epistemic and aleatoric uncertainties (§§ 3.1, 3.2). The result of the assimilation is the analysis state (§§ 3.4, 3.3), which provides a statistically optimal estimate of the physical quantity that we wish to predict, which is the truth ${\boldsymbol {\mathit {d}}}^{t }_k$. As explained in § 2.1, the truth is unknown thus we define the post-processed acoustic pressure as the presumed ground truth (figure 2). We drop the time subscript $k$ unless it is necessary for clarity.

3.1. Aleatoric and epistemic uncertainties

First, we discuss the statistical hypotheses on the aleatoric uncertainties. The aleatoric uncertainties contaminate the state and parameters as

(3.1a,b)\begin{equation} {\boldsymbol{\mathit{\phi}}} + {\boldsymbol{\mathit{\epsilon}}}_\phi = {\boldsymbol{\mathit{\phi}}}^{t}, \quad {\boldsymbol{\mathit{\alpha}}} + {\boldsymbol{\mathit{\epsilon}}}_\alpha = {\boldsymbol{\mathit{\alpha}}}^{t}, \end{equation}

where ${t}$ indicates the true quantity (which is unknown). The aleatoric uncertainties are modelled as Gaussian distributions

(3.2a,b)\begin{equation} {\boldsymbol{\mathit{\epsilon}}}_\phi \sim \mathcal{N}({\boldsymbol{\mathit{0}}}, \boldsymbol{\mathsf{C}}_{\phi\phi} ), \quad {\boldsymbol{\mathit{\epsilon}}}_\alpha \sim \mathcal{N}({\boldsymbol{\mathit{0}}}, \boldsymbol{\mathsf{C}}_{\alpha\alpha} ), \end{equation}

where $\mathcal {N}({\boldsymbol {\mathit {0}}}, \boldsymbol{\mathsf{C}})$ is a normal distribution with zero mean and covariance $\boldsymbol{\mathsf{C}}$. Second, we discuss model biases, which are epistemic uncertainties (Nóvoa et al. 2024). The model bias may arise from modelling assumptions and approximations in the operators $\mathcal {F}$ and $\mathcal {M}$. The true bias of a model is defined as

(3.3)\begin{equation} {\boldsymbol{\mathit{b}}}^{t} = {\boldsymbol{\mathit{d}}}^{t} -\mathbb{E}({\boldsymbol{\mathit{q}}}), \end{equation}

which is the expected difference between the true observable and the model observable. (If the model is unbiased, ${\boldsymbol {\mathit {d}}}^{t} =\mathbb {E}({\boldsymbol {\mathit {q}}})$.) The true bias is unknown (colloquially, an ‘unknown unknown’ (Nóvoa et al. 2024)) because we do not have access to the ground truth. Therefore, we need to estimate it from the information available from the data. We refer to ${\boldsymbol {\mathit {b}}}^{\dagger} = {\boldsymbol {\mathit {d}}}^{{\dagger}}-\mathbb {E}({\boldsymbol {\mathit {q}}})$ as the presumed true bias, and to the bias estimated from a model as ${\boldsymbol {\mathit {b}}}$. With this, the model equations, which define the first source of information on the system, are

(3.4) \begin{equation} \left.\begin{array}{c@{}} {\rm d}{\boldsymbol{\mathit{\phi}}}= \mathcal{F}(\boldsymbol{\mathit{\phi}}+{\boldsymbol{\mathit{\epsilon}}}_\phi, \boldsymbol{\alpha}+{\boldsymbol{\mathit{\epsilon}}}_\alpha ) {\rm d}t, \\ {}{\boldsymbol{\mathit{y}}} = \mathcal{M}({\boldsymbol{\mathit{\theta}}}, {\boldsymbol{\mathit{\phi}}}) + {\boldsymbol{\mathit{b}}}+ {\boldsymbol{\mathit{\epsilon}}}_q, \end{array} \right\} \end{equation}

where ${\boldsymbol {\mathit {y}}}$ is the bias-corrected model observable and ${\boldsymbol {\mathit {\epsilon }}}_q \sim \mathcal {N}({\boldsymbol {\mathit {0}}}, \boldsymbol{\mathsf{C}}_{qq})$. The set of equations (3.4) is not closed until we define an estimator for the model bias (§ 4).

The second source of information on the system is the data measured by the sensors. Experimental data are affected by both aleatoric noise and measurement shifts (§ 2.1). The true measurement shift is

(3.5)\begin{equation} {\boldsymbol{\mathit{b}}}_d^{t} = \mathbb{E}({\boldsymbol{\mathit{d}}}) - {\boldsymbol{\mathit{d}}}^{t}. \end{equation}

Because the ground truth is not known, the best estimate on the model bias is the presumed true measurement shift ${\boldsymbol {\mathit {b}}}_d^{\dagger} = {\boldsymbol {\mathit {d}}}^{\dagger} -\mathbb {E}({\boldsymbol {\mathit {d}}})$ and the estimated measurement shift is ${\boldsymbol {\mathit {b}}}_d$. With this, we define the measurements at a time instant as

(3.6)\begin{equation} {\boldsymbol{\mathit{d}}} + {\boldsymbol{\mathit{b}}}_d + {\boldsymbol{\mathit{\epsilon}}}_d = {\boldsymbol{\mathit{d}}}^{t}. \end{equation}

For clarity, we include a summary of the terminology in table 1. In the problem under investigation, ${\boldsymbol {\mathit {d}}}$ is the raw acoustic data, ${\boldsymbol {\mathit {d}}}^{{\dagger} }$ is the post-processed data and ${\boldsymbol {\mathit {b}}}_d$ is the non-zero mean of the raw data (see § 2.1). The aleatoric noise affects the measurement as ${\boldsymbol {\mathit {\epsilon }}}_d \sim \mathcal {N}({\boldsymbol {\mathit {0}}}, \boldsymbol{\mathsf{C}}_{dd})$. For simplicity, we assume that the measurement errors are statistically independent, i.e. $\boldsymbol{\mathsf{C}}_{dd}$ is a diagonal matrix with identical diagonal entries $\sigma _d$.

Table 1. Summary of the terminology.

Both the model bias and measurement shift are unknown a priori. To infer them, we analyse the residuals between the forecast and the observations, which are also known as innovations (Dee & Todling 2000; Haimberger 2007), which are defined as

(3.7)\begin{equation} {\boldsymbol{\mathit{i}}} = {\boldsymbol{\mathit{d}}} - {\boldsymbol{\mathit{q}}} \quad\Rightarrow\quad \mathbb{E}({\boldsymbol{\mathit{i}}}) = {\boldsymbol{\mathit{b}}}_d + {\boldsymbol{\mathit{b}}}. \end{equation}

Therefore, the expected innovation is the sum of the measurement and model biases as defined in (3.3) and (3.5). The relation between the innovations ${\boldsymbol {\mathit {i}}}$ and the biases ${\boldsymbol {\mathit {b}}}$ and ${\boldsymbol {\mathit {b}}}_d$ will be essential for the design of the bias estimator in § 4.

3.2. Augmented state-space formulation

We define the augmented state vector ${\boldsymbol {\mathit {\psi }}}=[{\boldsymbol {\mathit {\phi }}}; {\boldsymbol {\mathit {\alpha }}};{\boldsymbol {\mathit {q}}}]$, which comprises the state variables, the thermoacoustic parameters and the model observables to formally have a linear measurement operator $\boldsymbol{\mathsf{M}}$, which simplifies the derivation of DA methods (e.g. Nóvoa & Magri 2022). The augmented form of the model (3.4) yields

(3.8) \begin{equation} \left\{\begin{array}{@{}rcl} {\rm d}\begin{bmatrix} {\boldsymbol{\mathit{\phi}}}\\ {\boldsymbol{\mathit{\alpha}}}\\ {\boldsymbol{\mathit{q}}} \end{bmatrix} & = & \begin{bmatrix} \mathcal{F}({\boldsymbol{\mathit{\phi}}}+{\boldsymbol{\mathit{\epsilon}}}_\phi,{\boldsymbol{\mathit{\alpha}}}+{\boldsymbol{\mathit{\epsilon}}}_\alpha)\\ {\boldsymbol{\mathit{0}}}_{N_\alpha}\\ {\boldsymbol{\mathit{0}}}_{N_q} \end{bmatrix} {\rm d}t\\ {\boldsymbol{\mathit{y}}} & = & {\boldsymbol{\mathit{q}}} + {\boldsymbol{\mathit{b}}} + {\boldsymbol{\mathit{\epsilon}}}_{q} \end{array} \right. \;\leftrightarrow \; \left\{ \begin{array}{@{}rcl} {\rm d}{\boldsymbol{\mathit{\psi}}} & = & \boldsymbol{\mathsf{F}}\left({\boldsymbol{\mathit{\psi}}} +{\boldsymbol{\mathit{\epsilon}}}_\psi\right){\rm d}t \\ {\boldsymbol{\mathit{y}}} & = & \boldsymbol{\mathsf{M}}{\boldsymbol{\mathit{\psi}}} + {\boldsymbol{\mathit{b}}} + {\boldsymbol{\mathit{\epsilon}}}_{q} \end{array} \right. , \quad\forall\,t \neq {t_d}, \end{equation}

where $t_d$ are the times when the assimilation is performed, $\boldsymbol{\mathsf{F}}({\boldsymbol {\mathit {\psi }}})$ and ${\boldsymbol {\mathit {\epsilon }}}_\psi$ are the augmented nonlinear operator and aleatoric uncertainties, respectively; $\boldsymbol{\mathsf{M}} = [\boldsymbol{\mathsf{0}}~|~\mathbb {I}_{N_q}]$ is the linear measurement operator, which consists of the vertical concatenation of a matrix of zeros, $\boldsymbol{\mathsf{0}}\in \mathbb {R}^{N_q\times (N_\phi +N_\alpha )}$, and the identity matrix, $\mathbb {I}_{N_q}\in \mathbb {R}^{N_q\times N_q}$; and ${\boldsymbol {\mathit {0}}}_{N_\alpha }$ and ${\boldsymbol {\mathit {0}}}_{N_q}$ are vectors of zeros (because the parameters are constant in time, and ${\boldsymbol {\mathit {q}}}$ is not integrated in time but it is only computed at the analysis step).

3.3. Stochastic ensemble framework

Under the Gaussian assumption, the inverse problem of finding the states, ${\boldsymbol {\mathit {\phi }}}$, and parameters, ${\boldsymbol {\mathit {\alpha }}}$, given some observations, ${\boldsymbol {\mathit {d}}}$, would be solved by the Kalman filter equations if the operator $\mathcal {F}$ were linear (Kalman 1960). Thermoacoustic oscillations, however, have nonlinear dynamics (see (2.1)). Stochastic ensemble methods are suitable for nonlinear systems because they do not require to propagate the covariance, in contrast to other sequential methods, e.g. the extended Kalman filter (Evensen 2009; Nóvoa & Magri 2022). Stochastic ensemble filters track in time $m$ realizations of the augmented state ${\boldsymbol {\mathit {\psi }}}_j$ to estimate the mean and covariance, respectively

(3.9a)$$\begin{gather} \mathbb{E}({\boldsymbol{\mathit{\psi}}})\approx\bar{{\boldsymbol{\mathit{\psi}}}}=\dfrac{1}{m}\sum^m_{j=1}{{\boldsymbol{\mathit{\psi}}}_j} \end{gather}$$

(3.9b)$$\begin{gather}\boldsymbol{\mathsf{C}}_{\psi\psi} = \begin{bmatrix} \boldsymbol{\mathsf{C}}_{\phi\phi} & \boldsymbol{\mathsf{C}}_{\phi\alpha} & \boldsymbol{\mathsf{C}}_{\phi q} \\ \boldsymbol{\mathsf{C}}_{\alpha \phi} & \boldsymbol{\mathsf{C}}_{\alpha\alpha} & \boldsymbol{\mathsf{C}}_{\alpha q} \\ \boldsymbol{\mathsf{C}}_{q\phi} & \boldsymbol{\mathsf{C}}_{q\alpha} & \boldsymbol{\mathsf{C}}_{qq} \\ \end{bmatrix} \approx\dfrac{1}{m-1}\sum^m_{j=1}({\boldsymbol{\mathit{\psi}}}_j-\bar{{\boldsymbol{\mathit{\psi}}}})\otimes({\boldsymbol{\mathit{\psi}}}_j-\bar{{\boldsymbol{\mathit{\psi}}}}), \end{gather}$$

where $\otimes$ is the dyadic product. Because the forecast operator $\mathcal {F}$ is nonlinear, the Gaussian prior may not remain Gaussian after the model forecast, and $\mathbb {E}(\mathcal {F}({\boldsymbol {\mathit {\psi }}}))\neq \mathcal {F}(\bar {{\boldsymbol {\mathit {\psi }}}})$. However, the time between analyses $\Delta t_d$ is assumed small enough such that the Gaussian distribution is not significantly distorted (Evensen 2009; Yu, Juniper & Magri 2019b).

Finally, we approximate the model bias as ${\boldsymbol {\mathit {b}}}\approx {\boldsymbol {\mathit {d}}}^{t} - \boldsymbol{\mathsf{M}}\bar {{\boldsymbol {\mathit {\psi }}}}$; thus, the sum of the biases is approximately equal to the mean of the innovations, i.e.

(3.10)\begin{equation} \bar{{\boldsymbol{\mathit{i}}}} = {\boldsymbol{\mathit{d}}} - \boldsymbol{\mathsf{M}}\bar{{\boldsymbol{\mathit{\psi}}}} \approx {\boldsymbol{\mathit{b}}}_d + {\boldsymbol{\mathit{b}}}. \end{equation}

3.4. The r-EnKF

The objective function in a bias-regularized ensemble DA framework contains three norms (Nóvoa et al. 2024)

(3.11)\begin{equation} \mathcal{J}({\boldsymbol{\mathit{\psi}}}_j) = \left\|{\boldsymbol{\mathit{\psi}}}_j-{\boldsymbol{\mathit{\psi}}}_j^{f}\right\|^2_{\boldsymbol{\mathsf{C}}^{{f}^{{-}1}}_{\psi\psi}} + \left\|{{\boldsymbol{\mathit{y}}}}_j-{\boldsymbol{\mathit{d}}}_j\right\|^2_{\boldsymbol{\mathsf{C}}^{{-}1}_{dd}}+ \gamma\left\|{\boldsymbol{\mathit{b}}}_j\right\|^2_{\boldsymbol{\mathsf{C}}^{{-}1}_{bb}} \quad \mathrm{for}\ \ j=0,\dots,m-1, \end{equation}

where the superscript ‘f’ indicates ‘forecast’, the operator $\|{\cdot }\|^2_{\boldsymbol{\mathsf{C}}^{-1}}$ is the $L_2$-norm weighted by the semi-positive definite matrix ${\boldsymbol{\mathsf{C}}^{-1}}$, $\gamma \geq 0$ is a user-defined bias regularization factor and ${\boldsymbol {\mathit {b}}}_j$ is the model bias of each ensemble member. For simplicity, we define the bias in the ensemble mean, i.e. ${\boldsymbol {\mathit {b}}}_j = {\boldsymbol {\mathit {b}}}$ for all $j$ (which we estimate with an ESN in § 4) (Nóvoa et al. 2024). From left to right, the norms on the right-hand side of (3.11) measure (1) the spread of the ensemble prediction, (2) the distance between the bias-corrected estimate and the observables and (3) the model bias. The analytical solution of the r-EnKF in (3.12), which globally minimizes the cost function (3.11) with respect to ${\boldsymbol {\mathit {\psi }}}_j$, is

(3.12a) \begin{equation} {\boldsymbol{\mathit{\psi}}}_j^{a} = {\boldsymbol{\mathit{\psi}}}_j^{f} + \boldsymbol{\mathsf{K}}_r [(\mathbb{I}+ \boldsymbol{\mathsf{J}})^{\mathrm{T}}({\boldsymbol{\mathit{d}}}_j + {\boldsymbol{\mathit{b}}}_d - {\boldsymbol{\mathit{y}}}_j^{f}) - \gamma \boldsymbol{\mathsf{C}}_{dd}\boldsymbol{\mathsf{C}}^{{-}1}_{bb}\boldsymbol{\mathsf{J}}^{\mathrm{T}}{\boldsymbol{\mathit{b}}}], \quad j=0,\dots,m-1, \end{equation}

with

(3.12b) \begin{equation} \boldsymbol{\mathsf{K}}_r = \boldsymbol{\mathsf{C}}_{\psi\psi}^{f}\boldsymbol{\mathsf{M}}^{\mathrm{T}}[\boldsymbol{\mathsf{C}}_{dd} + (\mathbb{I}+ \boldsymbol{\mathsf{J}})^{\mathrm{T}}(\mathbb{I}+ \boldsymbol{\mathsf{J}})\boldsymbol{\mathsf{M}}\boldsymbol{\mathsf{C}}_{\psi\psi}^{f}\boldsymbol{\mathsf{M}}^{\mathrm{T}} + \gamma \boldsymbol{\mathsf{C}}_{dd}\boldsymbol{\mathsf{C}}^{{-}1}_{bb}\boldsymbol{\mathsf{J}}^{\mathrm{T}}\boldsymbol{\mathsf{J}}\boldsymbol{\mathsf{M}}\boldsymbol{\mathsf{C}}_{\psi\psi}^{f}\boldsymbol{\mathsf{M}}^{\mathrm{T}}]^{{-}1}, \end{equation}

where ‘a’ stands for ‘analysis’, i.e.  the optimal state of the assimilation; we assimilate each ensemble member with a different ${\boldsymbol {\mathit {d}}}_j\sim \mathcal {N}({\boldsymbol {\mathit {d}}}, \boldsymbol{\mathsf{C}}_{dd})$ to avoid covariance underestimation in ensemble filters (Burgers, Jan van Leeuwen & Evensen 1998); $\boldsymbol{\mathsf{K}}_r$ is the regularized Kalman gain matrix and $\boldsymbol{\mathsf{J}}= {\mathrm {d}} {\boldsymbol {\mathit {b}}}/ {\mathrm {d}}\boldsymbol{\mathsf{M}}{\boldsymbol {\mathit {\psi }}}$ is the Jacobian of the bias estimator. Formulae (3.12) generalize the analytical solutions of Nóvoa et al. (2024) to the case of a biased measurement. We prescribe $\boldsymbol{\mathsf{C}}_{dd} = \boldsymbol{\mathsf{C}}_{bb}$ because the model bias is defined in the observable space. We use $\gamma$ to tune the norm of the bias in (3.11) (Nóvoa et al. 2024). The optimal state and parameters are

(3.13) \begin{align} \begin{bmatrix} {\boldsymbol{\mathit{\phi}}}_j^{a}\\ {\boldsymbol{\mathit{\alpha}}}_j^{a} \end{bmatrix} &= \begin{bmatrix} {\boldsymbol{\mathit{\phi}}}_j^{f}\\ {\boldsymbol{\mathit{\alpha}}}_j^{f} \end{bmatrix} + \overbrace{ \begin{bmatrix} \boldsymbol{\mathsf{C}}_{\phi q}^{f}\\ \boldsymbol{\mathsf{C}}_{\alpha q}^{f} \end{bmatrix} \{\boldsymbol{\mathsf{C}}_{dd}+(\mathbb{I}+ \boldsymbol{\mathsf{J}})^{\mathrm{T}}(\mathbb{I}+ \boldsymbol{\mathsf{J}})\boldsymbol{\mathsf{C}}_{qq}^{f} +\gamma \boldsymbol{\mathsf{J}}^{\mathrm{T}}\boldsymbol{\mathsf{J}}\boldsymbol{\mathsf{C}}_{qq}^{f}\}^{{-}1} }^{\text{Regularized Kalman gain,}\, \boldsymbol{\mathsf{K}}_r} \ldots \nonumber\\ &\quad \dots\left[\left(\mathbb{I}+ \boldsymbol{\mathsf{J}}\right)^{\mathrm{T}}\overbrace{({\boldsymbol{\mathit{d}}}_j + {\boldsymbol{\mathit{b}}}_d - {\boldsymbol{\mathit{y}}}_j^{f})}^\text{Corrected innovation} - \gamma \boldsymbol{\mathsf{J}}^{\mathrm{T}}{\boldsymbol{\mathit{b}}}\right]. \end{align}

The r-EnKF defines a ‘good’ analysis from a biased model if the unbiased state ${\boldsymbol {\mathit {y}}}$ matches the truth, and the model bias ${\boldsymbol {\mathit {b}}}$ is small relative to the truth. The underlying assumptions of this work are that (i) our low-order model is qualitatively accurate such that the model bias ${\boldsymbol {\mathit {b}}}$ has a small norm, and (ii) the sensors are properly calibrated. In the limiting case when the assimilation framework is unbiased (i.e. ${\boldsymbol {\mathit {b}}}={\boldsymbol {\mathit {0}}}$ and ${\boldsymbol {\mathit {b}}}_d ={\boldsymbol {\mathit {0}}}$), the r-EnKF (3.13) becomes the bias-unregularized EnKF (Appendix A).

4. Reservoir computing for inferring model bias and measurement shift

To apply the r-EnKF (3.13), we must provide an estimate of the model bias and the measurement shift. We employ an ESN for this task. Echo state networks are suitable for real-time DA because (i) they are recurrent neural networks, i.e. they are designed to learn temporal dynamics in data; (ii) they are based on reservoir computing, hence, they are universal approximators (Grigoryeva & Ortega 2018); (iii) they are general nonlinear auto-regressive models (Aggarwal 2018); and (iv) training an ESN consists of solving a linear regression problem, which provides a global minimum without back propagation. In this work, we generalize the implementation of Nóvoa et al. (2024) to account for both model bias and measurement shift (§ 4.1). We propose a training dataset generation based on time-series correlation (§ 4.2).

4.1. Echo state network architecture and state-space formulation

We propose an ESN that simultaneously estimates ${\boldsymbol {\mathit {b}}}$ and ${\boldsymbol {\mathit {b}}}_d$. Importantly, ${\boldsymbol {\mathit {b}}}$ and ${\boldsymbol {\mathit {b}}}_d$ are not directly measurable in real time, but we have access to the innovation, which is linked to the biases through (3.10). Once we have information on the innovation and the model bias (the measurement shift), we can estimate the measurement shift (the model bias) through (3.10).

Figure 4(a) represents pictorially the proposed ESN at time $t_k$, with its three main components: (i) the input data, which are the mean analysis innovations, i.e. $\bar {{\boldsymbol {\mathit {i}}}}_k$; (ii) the reservoir, which is a high-dimensional state characterized by the sparse reservoir matrix $\boldsymbol{\mathsf{W}}\in \mathbb {R}^{N_{r}\times N_{r}}$ and vector ${\boldsymbol {\mathit {r}}}_k\in \mathbb {R}^{N_{r}}$ ($N_r\gg N_q$ is the number of neurons in the reservoir states); and (iii) the outputs, which are the mean innovation and the model bias at the next time step, i.e. $\bar {{\boldsymbol {\mathit {i}}}}_{k+1}$ and ${\boldsymbol {\mathit {b}}}_{k+1}$. The inputs to the network are a subset of the output (i.e. the innovation only). The sparse input matrix $\boldsymbol{\mathsf{W}}_{in}\in \mathbb {R}^{N_{r}\times (N_q+1)}$ maps the physical state into the reservoir, and the output matrix $\boldsymbol{\mathsf{W}}_{out}\in \mathbb {R}^{ N_q\times (N_r+1)}$ maps the reservoir state back to the physical state. Furthermore, figure 4(a) shows the two forecast settings of the ESN: open loop, which is performed when data on the innovations are available, and the closed loop, in which the ESN runs autonomously using the output as the input in the next time step. Figure 4(b) shows the unfolded architecture starting at time $t_k$ when observations of the innovations are available (i.e. at the analysis step). At this time, the reservoir is re-initialized with the analysis innovations such that the input to the ESN is $\bar {{\boldsymbol {\mathit {i}}}}^{a}_{k}={\boldsymbol {\mathit {d}}}_k-\boldsymbol{\mathsf{M}}\bar {{\boldsymbol {\mathit {\psi }}}}_k^{a}$. From this, the ESN estimates at the next time step $t_{k+1}$ the model bias and the innovation $\bar {{\boldsymbol {\mathit {i}}}}_{k+1}$, which is used as an initial condition in the subsequent forecast step. Mathematically, the equations that forecast the model bias and mean innovation in time are

(4.1)\begin{equation} \left.\begin{gathered} {}[{\boldsymbol{\mathit{b}}}_{k+1}; \bar{{\boldsymbol{\mathit{i}}}}_{k+1}] = \boldsymbol{\mathsf{W}}_{out}\left[{\boldsymbol{\mathit{r}}}_{k+1}; 1\right],\\ {\boldsymbol{\mathit{r}}}_{k+1} = \textrm{tanh}\left(\sigma_{in}\boldsymbol{\mathsf{W}}_{in}\left[\bar{{\boldsymbol{\mathit{i}}}}^\star_k\odot {\boldsymbol{\mathit{g}}}; \delta_r\right]+ \rho\boldsymbol{\mathsf{W}}{\boldsymbol{\mathit{r}}}_{k}\right), \end{gathered}\right\} \end{equation}

where $\bar {{\boldsymbol {\mathit {i}}}}^\star =\bar {{\boldsymbol {\mathit {i}}}}^{a}$ in open loop and $\bar {{\boldsymbol {\mathit {i}}}}^\star = \bar {{\boldsymbol {\mathit {i}}}}$ in closed loop; the ${\tanh ({\cdot })}$ operation is performed element-wise; the operator $\odot$ is the Hadamard product, i.e. component-wise multiplication; and ${\boldsymbol {\mathit {g}}}=[g_1;\dots ; g_{2N_q}]$ is the input normalizing term with $g_q^{-1} = \max {({\boldsymbol {\mathit {u}}}_q)}-\min {({\boldsymbol {\mathit {u}}}_q)}$, i.e. $g_q^{-1}$ is the range of the $q\mathrm {th}$ component of the training data $\boldsymbol{\mathsf{U}}$; $\delta _r$ is a constant used for breaking the symmetry of the ESN (we set $\delta _r=0.1$) (Huhn & Magri 2020). We define $\boldsymbol{\mathsf{W}}_{in}$ and $\boldsymbol{\mathsf{W}}$ as sparse and randomly generated, with a connectivity of 3 neurons (see Jaeger & Haas (2004) for details). We compute the matrix $\boldsymbol{\mathsf{W}}_{out}$ during the training.

Figure 4. Schematic representation of the proposed ESN architecture for model bias and innovation prediction. (a) Compact architecture showing the open-loop and closed-loop configurations; and (b) unfolded architecture starting at an analysis step is at $t_k$, in which there is one open-loop step followed by a closed-loop forecast. In open loop the input to the ESN is the analysis mean innovation $\bar {{\boldsymbol {\mathit {i}}}}^{a}_{k}={\boldsymbol {\mathit {d}}}_k-\boldsymbol{\mathsf{M}}\bar {{\boldsymbol {\mathit {\psi }}}}_k^{a}$, i.e. the difference between the raw acoustic pressure data and the analysis model estimate. The outputs from the ESN are (i) the innovation at the next time step $\bar {{\boldsymbol {\mathit {i}}}}_{k+1}$ (which becomes the input to the next time step in the closed-loop configuration); and (ii) the model bias ${\boldsymbol {\mathit {b}}}_{k+1}$, i.e. an estimate of the difference between the presumed acoustic state ${\boldsymbol {\mathit {d}}}_{k+1}^{\dagger}$ and the model estimate $\boldsymbol{\mathsf{M}}\bar{\boldsymbol{\psi}}_{k+1} $.

Lastly, the r-EnKF equations (3.12) require the definition of the Jacobian of the bias estimator at the analysis step. Because at the analysis step the ESN inputs the analysis innovations with an open-loop step (see figure 4b), the Jacobian of the bias estimator is equivalent to the negative Jacobian of the ESN in open-loop configuration (see Appendix B).

4.2. Training the network

During training, the ESN is in an open-loop configuration (figure 4a). The inputs to the reservoir are the training dataset $\boldsymbol{\mathsf{U}}=[{\boldsymbol {\mathit {u}}}_0\,|\dots |\,{\boldsymbol {\mathit {u}}}_{N_{tr}-1}]$, in which each time component ${\boldsymbol {\mathit {u}}}_k = [{\boldsymbol {\mathit {b}}}_{u}(t_k); {{\boldsymbol {\mathit {i}}}}_{u}(t_k)]$, the subscript ‘$u$’ indicates training data, and the operator $[~|~]$ indicates horizontal concatenation. Although we have information from the experimental data to train the network, the optimal parameters of the thermoacoustic system are unknown. Thus, we do not know a priori the model bias and measurement shift. Selecting an appropriate training dataset is key to obtaining a robust ESN, which can estimate the model bias and innovations. We create a set of $L$ guesses on the bias and innovations from a single realization of the experimental data, which means that the ESN is not trained with the ‘true’ bias (which is an unknown quantity in real time). The training data are generated from the experimental data as detailed in Appendix C. In summary, the procedure is as follows.

1. (i) Take measurements for a training time window $t_{tr}$ of acoustic pressure data $\boldsymbol{\mathsf{D}}_{u}$ and estimate $\boldsymbol{\mathsf{D}}^{{\dagger} }_{u}$ by applying a Butterworth filter to the raw data $\boldsymbol{\mathsf{D}}_{u}$ (see § 2).

2. (ii) Generate $L$ model estimates $\boldsymbol{\mathsf{Q}}_{{u}, l}$ of the acoustic pressure from (2.4). Each $\boldsymbol{\mathsf{Q}}_{{u}, l}$ has a different set of parameters, which are uniformly randomly generated (the parameters’ ranges are reported in Appendix D).

3. (iii) Correlate in time each $\boldsymbol{\mathsf{Q}}_{{u}, l}$ with $\boldsymbol{\mathsf{D}}_{u}$.

4. (iv) Compute the model bias and innovations as $\boldsymbol{\mathsf{B}}_{{u},l}= \boldsymbol{\mathsf{D}}^{{\dagger} }_{u} - \boldsymbol{\mathsf{Q}}_{{u},l}$ and $\boldsymbol{\mathsf{I}}_{{u},l} = \boldsymbol{\mathsf{D}}_{u} - \boldsymbol{\mathsf{Q}}_{{u},l}.$

Finally, we apply data augmentation to improve the network's adaptivity in a real-time assimilation framework. The total number of training time series in the proposed training method is $2L$. Training the ESN consists of finding the elements in $\boldsymbol{\mathsf{W}}_{out}$, which minimize the distance between the outputs obtained from an open-loop forecast step of the training data (the features) to the training data at the following time step (the labels). This minimization is solved by ridge regression of the linear system (Lukoševičius 2012)

(4.2)\begin{equation} \left(\sum_{l = 0}^{2L-1}\boldsymbol{\mathsf{R}}^{\,}_l\boldsymbol{\mathsf{R}}_l^{\mathrm{T}} + \lambda\, \mathbb{I}_{N_r+1}\right)\boldsymbol{\mathsf{W}}_{out}^{\mathrm{T}} = \sum_{l=0}^{2L-1}\boldsymbol{\mathsf{R}}_l^{\,}\boldsymbol{\mathsf{U}}_l^{\mathrm{T}}, \end{equation}

where $\lambda$ is the Tikhonov regularization parameter and $\boldsymbol{\mathsf{R}} = [[{\boldsymbol {\mathit {r}}}_0; 1]\,|\dots |\,[{\boldsymbol {\mathit {r}}}_{N_{tr}-1}; 1]]$, with ${\boldsymbol {\mathit {r}}}_0={\boldsymbol {\mathit {0}}}$ and ${\boldsymbol {\mathit {r}}}_k$ are obtained with (4.1) using the innovations of the training set as inputs. The summations over $2L$ can be performed in parallel to minimize computational costs. Finally, we tune the hyperparameters of the ESN, i.e. the spectral radius $\rho$, the input scaling $\sigma _{in}$ and the Tikhonov regularization parameter $\lambda$. We use a recycle validation strategy with Bayesian optimization for the hyperparameter selection (see Racca & Magri 2021).

5. Implementation

To summarize, the following four stages are necessary to design the real-time digital twin.

(1) Initialization: the ensemble ${\boldsymbol {\mathit {\psi }}}_j\ \mathrm {for}\ j = 1,\ldots,m$ and the ESN are initialized using the parameters reported in Appendix D. (2) Forecast: time march in parallel each ensemble member (3.8), i.e. the system of thermoacoustic equations and parameters, and the ESN according to (4.1) until observation data become available. (3) Analysis: apply the r-EnKF (3.12), which obtains the optimal combination between the unbiased model estimate and the observation data. (4) Re-initialization: update the state and parameters of the ensemble with the analyses (i.e. ${\boldsymbol {\mathit {\phi }}}^{a}$ and ${\boldsymbol {\mathit {\alpha }}}^{a}$), and the ESN with the analysis innovation (i.e. $\bar {{\boldsymbol {\mathit {i}}}}^{a} = {\boldsymbol {\mathit {d}}} - \boldsymbol{\mathsf{M}}\bar {{\boldsymbol {\mathit {\psi }}}}^{a}$).

Steps (2)–(4) are repeated sequentially as data become available. Once the assimilation process has ended, we forecast further the ensemble and the ESN to analyse the extrapolation and generalization capability of the estimated state, parameters and biases.

6. Performance of the real-time digital twin

The real-time digital twin has four components: (i) data from sensors (§ 2.1), (ii) the physics-based low-order model (§ 2.2), (iii) the ESN to infer the model bias and measurement shift (§ 4), and (iv) the real-time and bias-regularized DA method (§ 3.4) to twin sensors’ data with the low-order model to infer the state and model parameters. We investigate the capability of the proposed real-time digital twin to predict autonomously the dynamics of azimuthal thermoacoustic oscillations. We analyse the digital twin's performance with respect to established methods: the bias-unregularized EnKF (Nóvoa & Magri 2022), which is suitable for real-time assimilation, and Langevin-based regression (Indlekofer et al. 2022), which is suitable for offline assimilation. The hyperparameters used to train the ESN and the assimilation parameters are reported in Appendix D.

6.1. Time-accurate prediction

We focus our analysis on the equivalence ratio $\varPhi =0.5125$, which is a thermoacoustically unstable system. Figure 5 shows the time evolution of the acoustic pressure at $\theta =60^\circ$. Prior to the assimilation, we train the ESN with training data over $t_{tr}=0.167$ s. We start the assimilation at $0.5\ {\rm s}$ to avoid using training data during the assimilation and to have the ensemble in the statistically stationary regime. In parallel, we initialize the ESN before the assimilation begins with 10 data points (i.e. the washout). After this, observations from the raw experimental data (red dots) are assimilated every $\Delta t_d=5.86\times 10^{-4}$ s, which corresponds to approximately $0.64$ data points per acoustic period. (The assimilation frequency was selected by down-sampling the raw data, which is recorded at 51.2 kHz, such that $\Delta t_d$ fulfils to the Shannon–Nyquist criterion.) At a time instant, the measurement (the data point, red dot) is assimilated into the model, which adaptively updates itself through the r-EnKF and bias estimators (§§ 3.4, 4). After that, the state and model parameters have been updated, the low-order model (2.5) evolves autonomously until the next data point becomes available. This process mimics a stream of data coming from sensors on the fly. We assimilate 600 measurements during the assimilation window of 0.35 s, the model runs autonomously without seeing more data ($t>0.85$ s).

Figure 5. Real-time digital twin at $\varPhi =0.5125$. Time evolution of the thermoacoustic pressure (in Pa) at $\theta =60^\circ$ using (a) the bias-unregularized EnKF and (b) the proposed r-EnKF. Comparison between the presumed ground truth (thick grey line), the raw data (thick red line), the prediction from the ensemble filter (cyan lines) and in (b) the bias-corrected mean estimate (navy dashed line). The close ups show the start and end of the assimilation window, which is indicated by the vertical dashed lines. The observations are show in red circles.

First, we analyse the performance of the bias-unregularized filter (EnKF) on state estimation (figure 5a) with the corresponding parameter inference shown in figure 6 (dashed lines). Although the system has a self-sustained oscillation, the bias-unregularized filter converges towards an incorrect solution, that is, a fixed point. On the other hand, the bias-regularized filter (r-EnKF) successfully learns the thermoacoustic model. Second, the r-EnKF filters out aleatoric noise and turbulent fluctuations. These fluctuations are the difference between the amplitudes of the raw data, which are the input to the digital twin, and the post-processed data, which are the hidden state that we wish to uncover and predict (the presumed acoustic state), of up to approximately 200 Pa (close-ups at the end of assimilation in figure 5). Third, after assimilation, the r-EnKF predicts the thermoacoustic limit cycle beyond the assimilation window. This means that the digital twin has learnt an accurate model of the system, which generalizes beyond the assimilation. The ensemble model observables ${\boldsymbol {\mathit {q}}}_j$ and bias-corrected ensemble mean $\bar {{\boldsymbol {\mathit {y}}}} = \bar {{\boldsymbol {\mathit {q}}}} + {\boldsymbol {\mathit {b}}}$ are almost identical at convergence in figure 5(b), which means that the magnitude of the model bias is small (approximately 20 Pa in amplitude; see figure 8). Fourth, the r-EnKF infers the optimal system's thermoacoustic parameters, which are seven in this model (figure 6). The ensemble parameters at $t=0$ in figure 6 are initialized on physical ranges from the physical parameters used in the literature, which are computed by Langevin regression in Indlekofer et al. (2022). From inspection of the optimal parameters, we can draw physical conclusions: (i) the parameters of the literature are not the optimal parameters of the model; (ii) the linear growth rate, $\nu$, angular frequency, $\omega$, the heat release strength weighted by the symmetry intensity, $c_2\beta$, and the phase of the reactive symmetry, $\varTheta _\epsilon$, do not significantly vary at regime, which means that these parameters are constants and do not depend on the state; and (iii) the nonlinear parameter, $\kappa$, the amplitude of the reactive symmetry, $\epsilon$, and the direction of asymmetry, $\varTheta _\beta$ have temporal variations that follow the modulation of the pressure signal envelope (figure 6). Physically, this means that the optimal deterministic system that represents azimuthal instabilities is a time-varying parameter system (Durbin & Koopman 2012). By inferring time-varying parameters, we derive a deterministic system that does not require stochastic modelling. We further analyse the thermoacoustic parameters for all the equivalence ratios in figures 10 and 11.

Figure 6. Thermoacoustic parameters with the bias-unregularized EnKF (dashed) and the r-EnKF (solid). The lines and shaded areas indicate the ensemble mean and standard deviation, respectively. Here $\varPhi =0.5125$.

To summarize, the presence of model bias and measurement shift makes the established bias-unregularized EnKF fail to uncover and predict the physical state from noisy data. In contrast, the bias-regularized filter (r-EnKF) infers the state, parameters, model bias and measurement shift, which enable a time accurate and physical prediction beyond the assimilation window.

6.2. Statistics and biases

We analyse the uncertainty and the statistics of the time series (§ 6.1) generated by the real-time digital twin. To do so, we define the normalized r.m.s. error of two time series ${\boldsymbol {\mathit {w}}}, {\boldsymbol {\mathit {z}}}$ with $N_q$ dimensions and $N_t$ time steps as

(6.1)\begin{equation} \mathrm{R.m.s.}({\boldsymbol{\mathit{w}}}, {\boldsymbol{\mathit{z}}}) = \sqrt{\dfrac{\sum_{q}\sum_{k} \left({w}_q(t_k)- {z}_q(t_k)\right)^2}{\sum_{q} \sum_{k}\left({w}_q(t_k)\right)^2}}, \quad \mathrm{for} \ \begin{array}{l} q=0,\dots, N_q-1, \\ k=0,\dots,N_t-1. \end{array} \end{equation}

Figure 7 shows the r.m.s. errors at four critical time instants in the assimilation process: before the data are assimilated (first column), which corresponds to the initialization of the digital twin; at the start and end of the DA window (second and third columns, respectively); and after the DA (fourth column), i.e. when the digital twin evolves autonomously to predict unseen dynamics (generalization). The inference of the model bias and measurement shift is key to obtaining a small generalization error. Although the bias-unregularized EnKF converges to an unphysical solution with a large r.m.s., the bias-regularized filter converges to a physical solution with a small r.m.s. As the assimilation progresses, the ESN improves the prediction of the model bias because the r.m.s. of the bias-corrected prediction (figure 7b, cyan) is smaller than the model estimate (figure 7b, navy). This is further evidenced by analysing the histograms of the time series for the four azimuthal locations (figure 8). The digital twin (bias-corrected solution, navy histograms) converges to the expected distribution of the acoustic pressure (presumed truth) (black histograms) despite the fact that the assimilated data are substantially contaminated by noise and turbulent fluctuations (red histograms). The bias-corrected histograms have a zero mean, which means that the ESN in the digital twin has correctly inferred both the model bias and the measurement shift. Both biases are shown in figure 9. The measurement shift, whose presumed true value is known a priori ${\boldsymbol {\mathit {b}}}_d^{\dagger} = \langle {\boldsymbol {\mathit {d}}}-{\boldsymbol {\mathit {d}}}^{\dagger} \rangle$ (where the brackets $\langle ~\rangle$ indicate time average), is exactly inferred by the ESN. On the other hand, the presumed true value of the model bias is not known a priori, but it is presumed as ${\boldsymbol {\mathit {b}}}^{\dagger} = {\boldsymbol {\mathit {d}}}^{\dagger} - {\boldsymbol {\mathit {q}}}^{\dagger}$, where ${\boldsymbol {\mathit {q}}}^{\dagger}$ is the presumed true model estimate, which is not known a priori. The digital twin improves the knowledge that we have on the model bias by correcting the presumed model bias.

Figure 7. Normalized r.m.s. errors between the presumed ground truth and the model prediction (cyan) and between the presumed ground truth and the bias-corrected prediction (navy) with (a) the bias-unregularized EnKF and (b) the r-EnKF. The filled histograms show the r.m.s. of each ensemble member, the thick lines show the r.m.s. of the ensemble mean and the vertical dashed line indicates the mean of the r.m.s. error. The errors are computed at four stages of the assimilation: $t\in [0.49, 0.50]$ (before DA), $t\in [0.50, 0.51]$ (start of DA), $t\in [0.84, 0.85]$ (end of DA) and $t\in [0.85, 0.86]$ (after DA). Here $\varPhi =0.5125$.

Figure 8. Histograms of the acoustic pressure after assimilation with (a) the bias-unregularized EnKF and (b) the r-EnKF, at the four observed azimuthal locations. Comparison between the presumed ground truth (black), the observations from raw data (red), the ensemble prediction (filled cyan) and its mean (dashed teal), and bias-corrected ensemble predictions (filled navy) and its mean (dashed light blue). The vertical lines indicate the mean of each of the distributions. Here $\varPhi =0.5125$.

Figure 9. Model bias and measurement shift estimates after assimilation with the r-EnKF. Comparison between the ESN inference of (a) the model bias (dashed black) and the presumed model bias (thick grey); and (b) the measurement shift estimate (dashed black), difference between the raw and post-processed data (thin grey), and its time average, $\langle {\cdot }\rangle$ (thick grey). Here $\varPhi =0.5125$.

To summarize, the r-EnKF generalizes to unseen scenarios and extrapolates correctly in time. Key to the robust performance is the inference of the biases in the model and measurements.

6.3. Physical parameters and comparison with the literature

We analyse the system's physical parameters (§ 2.2 and (2.4)) for all the available equivalence ratios. We train an ESN for each equivalence ratio. The DA parameters (assimilation window, frequency, ensemble size, etc.) are the same in all tests (Appendix D). We compare the low-order model parameters (i) inferred by the bias-unregularized EnKF, (ii) inferred by the r-EnKF, and (iii) obtained with the state-of-the-art Langevin regression by Indlekofer et al. (2022). The linear growth rate, $\nu$, and the heat source strength weighted by the symmetry intensity, $c_2\beta$ are shown in figure 10; and the remaining parameters are shown in figure 11. The specific value parameters are listed in Appendix E.

Figure 10. Comparison of the estimated values of the two stability parameters $c_2\beta$ (circles) and $\nu$ (triangles) after DA, obtained with the bias-unregularized EnKF (black) and the r-EnKF with different bias regularization parameters (colour maps). Lighter colours indicate stronger regularization; larger marker sizes indicate more accurate solutions, i.e. smaller r.m.s. errors; and the error bars represent the ensemble standard deviation. The dashed lines corresponds to the parameters identified by offline Langevin regression (Indlekofer et al. 2022).

Figure 11. Same as figure 10 for the remaining physical parameters.

First, the r-EnKF infers physical parameters with a small generalization error (the r.m.s. errors are small, depicted with large markers). Large values of the bias regularization, i.e. $\gamma \gtrsim 1$, provide the most accurate set of parameters. This physically means that the low-order model is a qualitatively accurate model because the bias has a small norm (which corresponds to a large $\gamma$). The bias-unregularized EnKF is outperformed in all cases. (The case with $\varPhi =0.4875$ is a stable configuration with a trivial fixed point, which is why the parameters have a larger r.m.s.) Second, we analyse the uncertainty of the parameters, which is shown by an error bar. The height of the error bar is the ensemble standard deviation. The linear growth rate, $\nu$, and the angular frequency, $\omega$, are inferred with small uncertainties. Physically, this means that the model predictions are markedly sensitive to small changes in the linear growth rate and angular frequency. The prediction's sensitivity to the parameters change with the equivalence ratio. Physically, the large sensitivity of the linear growth rate (i.e. the system's linear stability) to the operating condition is a characteristic feature of thermoacoustic instabilities, in particular, in annular combustors, in which eigenvalues are degenerate (e.g. Magri et al. 2023). Third, there exists multiple combinations of parameters that provide a similar accuracy. For example, in the case $\varPhi =0.5125$, there are different combinations of $\varTheta _\beta, \epsilon, c_2\beta$, which yield a small r.m.s. Fourth, we compare the digital twin with the offline Langevin regression of Indlekofer et al. (2022) (horizontal dashed lines in figures 10 and 11). The inferred parameters are physically meaningful because they lie within value ranges that are similar to the parameters of the literature. The digital twin simultaneously infers all the parameters in real time, which overcomes the limitations of Langevin regression, which needs to be performed on one parameter at a time and offline. The parameters presented in figures 10 and 11 depend linearly on $\varPhi$. These linear dependencies are approximations of the ground truth and were inferred from the identified set of parameters at each $\varPhi$ (see figure 9 in Indlekofer et al. 2022). When tuned for a specific $\varPhi$, the offline parameters lead to state statistics that more closely match the observables. The digital twin parameters are optimal because they are the minimizers of (3.11). The long-term statistics of the case $\varPhi =0.5625$, which lives on a generalized Bloch sphere (Magri et al. 2023), is further analysed with the quaternion ansätz of Ghirardo & Bothien (2018). Figure 12 shows the acoustic state for $\varPhi =0.5625$ from the raw data (panel b), and the proposed digital twin (panel c). The long-term statistics are obtained by forecasting the model (without assimilation) using the identified parameters (reported in Appendix E). The real-time digital twin infers a set of physical parameters, which correctly capture the physical state of the azimuthal acoustic mode. The nature angle has a bimodal distribution with $|\chi |<{\rm \pi} /4$, which means that the thermoacoustic instability is a mixed mode.

Figure 12. Acoustic state for $\varPhi =0.5625$. (a) Bloch sphere for quaternion representation of the acoustic pressure $p(\theta, t) = A \cos {(\vartheta - \theta )} \cos {\chi } \cos {(\omega t + \varphi )} + A \sin {(\vartheta - \theta )} \sin {\chi } \sin {(\omega t + \varphi )}$, where $A$ is the slow-varying amplitude, i.e. the envelope of the acoustic pressure, $\varphi$ is the slow temporal phase drift, $\vartheta$ is the position of the anti-nodal line and $\chi$ is the nature angle (Ghirardo & Bothien 2018; Magri et al. 2023). (b) State from the raw data. (c) State predicted by the proposed real-time digital twin.

6.4. Generalizability

In this section we propose a method for generalizing the digital twin to unseen experimental data, e.g. a different equivalence ratio. In §§ 6.1–6.3 we train a different ESN for each equivalence ratio. Here, we test the real-time digital twin framework using an unified ESN, i.e. we train the ESN with data for $\varPhi ={0.4875, 0.5125, 0.5375}$ but not for the $\varPhi _{test}=0.5625$ case. The training data generation is identical to that detailed in Appendix C, but we combine the training data for $\varPhi ={0.4875, 0.5125, 0.5375}$ to train one ESN, which has the same characteristics as those used in the previous analyses (Appendix D). We increased the ensemble size to $m=40$ to account for the larger uncertainty and variability in the dynamics. Figure 13 shows the inferred linear growth rate $\nu$ and the heat source strength $c_2\beta$ for all $\varPhi$. Notably, the digital twin for the equivalence ratio $\varPhi _{test}=0.5625$, which was unseen by the ESN, converges to similar parameters to those in figure 10. We further analyse the test case $\varPhi _{test}=0.5625$ in figure 14, which shows the histograms of the acoustic pressure after assimilation. The ESN infers the correct measurement shift to give a zero mean prediction of the pressure, and a model bias that reduces the distance between the ensemble state estimate and the presumed truth. Overall, we conclude that the unified ESN successfully generalizes to data with unseen dynamics.

Figure 13. Generalizability study. Estimated values of the two stability parameters $c_2\beta$ (circles) and $\nu$ (triangles) after DA with the r-EnKF using the same ESN, which is not trained with data for $\varPhi =0.5625$. Lighter colours indicate a stronger bias regularization parameter; larger marker sizes indicate more accurate solutions, i.e. smaller r.m.s. errors; and the error bars represent the ensemble standard deviation. The dashed lines correspond to the parameters identified by offline Langevin regression (Indlekofer et al. 2022).

Figure 14. Generalizability study. Histograms of the acoustic pressure at $\varPhi =0.5625$ after assimilation with the r-EnKF, at the four observed azimuthal locations. The ESN used in the digital twin is not trained with data for $\varPhi =0.5625$. Comparison between the presumed ground truth (black), the observations from raw data (red), the ensemble prediction (filled cyan) and its mean (dashed teal), and bias-corrected ensemble predictions (filled navy) and its mean (dashed light blue). The vertical lines indicate the mean of each of the distributions.

7. Conclusions

We develop a real-time digital twin of azimuthal thermoacoustic oscillations of a hydrogen-based combustor for uncovering the physical state from raw data, and predicting the nonlinear dynamics in time. We identify the four ingredients required to design a robust real-time digital twin as follows.

1. (i) In a real-time context, we need to work with raw data as they come from sensors, therefore, we do not pre-process the data. The data we utilize is raw (the data contain both environmental and instrumental noise, and turbulent fluctuations) and sparse (the data are collected by four microphones).

2. (ii) A low-order model, which is deterministic, qualitatively accurate, computationally cheap and physical.

3. (iii) An estimator of both the model bias, which originates from the modelling assumptions and approximation of a low-order model, and measurement shift, which originates from microphones recording the total pressure (which has a non-zero mean) instead of the acoustic pressure (which has zero mean).

4. (iv) A statistical DA method that optimally combines the two sources of information on the system (i)–(ii) to improve the prediction on the physical states by updating the physical parameters every time that sensors’ data become available (on the fly or in real time).

We propose a real-time DA framework to infer the acoustic pressure, physical parameters, model bias and measurements shift simultaneously, which is the r-EnKF. We find an analytical solution that minimizes the bias-regularized DA cost function (3.11). We propose a reservoir computer (an ESN) and a training strategy to infer both the model bias and measurement shift, without making assumptions on their functional forms a priori. The real-time digital twin is applied to a laboratory hydrogen-based annular combustor for a variety of equivalence ratios. The data are treated as if they came from sensors on the fly, i.e. the pressure measurements are assimilated at a time step and then disregarded until the next pressure measurement becomes available. First, after assimilation of raw data (at 1707 Hz, i.e. approximately 0.6 data points per acoustic period), the model learns the correct physical state and optimal parameters. The real-time digital twin autonomously predicts azimuthal dynamics beyond the assimilation window, i.e. without seeing more data. This is in stark contrast with state-of-the-art methods based on the bias-unregularized methods, which do not perform well in the presence of model bias and/or measurement shift. Second, the digital twin uncovers the physical acoustic pressure from the raw data. Because the physical mechanisms are constrained in the low-order model, which originates from conservation laws, the digital twin acts as a physics-based filter, which removes aleatoric noise and turbulent fluctuations. Third, physically we find that azimuthal oscillations are governed by a time-varying parameter system, which generalizes existing models that have constant parameters and capture only slow-varying variables. Fourth, we find that the key parameters that influence the dynamics are the linear growth rate and the angular frequency. Furthermore, the digital twin generalizes to equivalence ratios at which current modelling approaches are not accurate. Fifth, to generalize the digital twin framework, we train an ESN with data from three of the four available equivalence ratios and test the real-time digital twin performance in the unseen scenario. This unified ESN successfully estimates the model bias and measurement shift in unseen thermoacoustic dynamics with the regularized bias-aware DA framework. This work opens new opportunities for real-time digital twinning of multi-physics problems.