2.1. Experimental set-up and data

We use the experimental data of Faure-Beaulieu et al. (Reference Faure-Beaulieu, Indlekofer, Dawson and Noiray2021a) and Indlekofer et al. (Reference Indlekofer, Faure-Beaulieu, Dawson and Noiray2022). 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. Reference Indlekofer, Faure-Beaulieu, Dawson and Noiray2022). 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. (Reference Faure-Beaulieu, Indlekofer, Dawson and Noiray2021b).

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 Reference Faure-Beaulieu and Noiray2020; Indlekofer et al. Reference Indlekofer, Faure-Beaulieu, Dawson and Noiray2022)

(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. Reference Indlekofer, Faure-Beaulieu, Dawson and Noiray2022); $\dot {q}$ is the coherent component of the fluctuations of the heat release rate (Noiray Reference Noiray2017), 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 Reference Lieuwen2003; Noiray Reference Noiray2017). 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. Reference Noiray, Bothien and Schuermans2011; Noiray & Schuermans Reference Noiray and Schuermans2013b; Magri et al. Reference Magri, Schmid and Moeck2023)

(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. Reference Indlekofer, Faure-Beaulieu, Dawson and Noiray2022). 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.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. Reference Nóvoa, Racca and Magri2024). 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. Reference Nóvoa, Racca and Magri2024)) 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 Reference Dee and Todling2000; Haimberger Reference Haimberger2007), 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 Reference Nóvoa and Magri2022). 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 Reference Kalman1960). 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 Reference Evensen2009; Nóvoa & Magri Reference Nóvoa and Magri2022). 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 Reference Evensen2009; Yu, Juniper & Magri Reference Yu, Juniper and Magri2019b).

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. Reference Nóvoa, Racca and Magri2024)

(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. Reference Nóvoa, Racca and Magri2024). 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 Reference Burgers, Jan van Leeuwen and Evensen1998); $\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. (Reference Nóvoa, Racca and Magri2024) 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. Reference Nóvoa, Racca and Magri2024). 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.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 Reference Huhn and Magri2020). We define $\boldsymbol{\mathsf{W}}_{in}$ and $\boldsymbol{\mathsf{W}}$ as sparse and randomly generated, with a connectivity of 3 neurons (see Jaeger & Haas (Reference Jaeger and Haas2004) 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 Reference Lukoševičius2012)

(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 Reference Racca and Magri2021).

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. (Reference Indlekofer, Faure-Beaulieu, Dawson and Noiray2022). 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 Reference Durbin and Koopman2012). 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. (Reference Indlekofer, Faure-Beaulieu, Dawson and Noiray2022). 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. Reference Indlekofer, Faure-Beaulieu, Dawson and Noiray2022).

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. Reference Magri, Schmid and Moeck2023). 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. (Reference Indlekofer, Faure-Beaulieu, Dawson and Noiray2022) (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. Reference Indlekofer, Faure-Beaulieu, Dawson and Noiray2022). 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. Reference Magri, Schmid and Moeck2023), is further analysed with the quaternion ansätz of Ghirardo & Bothien (Reference Ghirardo and Bothien2018). 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 Reference Ghirardo and Bothien2018; Magri et al. Reference Magri, Schmid and Moeck2023). (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. Reference Indlekofer, Faure-Beaulieu, Dawson and Noiray2022).

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.