Bias-aware thermoacoustic data assimilation

Andrea Nóvoa University of Cambridge Cambridge University Engineering Department, Trumpington St, Cambridge CB2 1PZ, UK

Alberto Racca University of Cambridge Cambridge University Engineering Department, Trumpington St, Cambridge CB2 1PZ, UK

Luca Magri 1

Imperial College London, Aeronautics Department, Exhibition Road, London SW7 2AZ, UK Cambridge University Engineering Department, Trumpington St, Cambridge CB2 1PZ, UK The Alan Turing Institute, 96 Euston Rd, London NW1 2DB, UK Institute for Advanced Study, TU Munich, Lichtenbergstraße 2a 85748 Garching, Germany (visiting)

ABSTRACT The occurrence and amplitude of nonlinear thermoacoustic instabilities can be quickly estimated with low-order models. Low-order models, however, contain model errors, i.e., the equations do not capture all the physical mechanisms. From a statistical inference point of view, we say that low-order models are biased. We propose a data-assimilation methodology that can simultaneously infer the acoustic state, model parameters and model error from reference data. We propose reservoir computing (echo state networks) to represent the model bias. The echo state network is combined with an ensemble square-root Kalman filter to perform, in real time, (i) the estimation of the acoustic pressure and velocity, (ii) the inference of the flame parameters, and (iii) the learning of the model bias. The proposed methodology is tested on a time-delayed low-order model, with synthetic experimental data from a higher-order model with a flame. This work opens up new possibilities for assimilating data for cheap low-order models to self-correct on-the-fly.

1. INTRODUCTION

Thermoacoustic oscillations arise in gas turbines or aeroengines due to the constructive interaction between a heat source (e.g. a flame), acoustics and hydrodynamics [1]. When the acoustic pressure and the flame are su ffi ciently in phase, the acoustic modes grow in amplitude, which gives rise to self-sustained oscillations [2]. These can produce noise, vibrations and may be critical to the turbomachinery [3]. The dynamics of thermoacoustics are nonlinear and extremely sensitive to perturbations to the operating conditions [4, 5], which hampers the manufacturing and control of thermoacoustically stable devices. The development of optimal design and control strategies relies on the mathematical model used to represent the dominant physical mechanisms [6]. Thermoacoustics

1 l.magri@imperial.ac.uk

a slaty. inter.noise 21-24 AUGUST SCOTTISH EVENT CAMPUS O ¥, ? GLASGOW

can be modelled with qualitative accurate low-order models. In contrast to high-fidelity simulations, low-order models are computationally inexpensive, however, the downside is the lack of quantitatively accuracy due to the simplifying assumptions. We develop a data assimilation method that improves the accuracy of thermoacoustic low–order models on the fly. Traverso & Magri [7] introduced data assimilation in thermoacoustics through a variational approach. Variational methods require batches of data for a minimisation problem, making them not suitable for on-the-fly inference [8]. Nonetheless, as accurate thermoacoustic experimental data becomes more accessible in real time [9], there is a growing interest in real- time inference methods, such as statistical sequential data assimilation, which bypass the need for post-processing and data storage. These methods correct the model predictions when reference data becomes available by performing a Bayesian update, which provides the maximum a posteriori estimate of the prediction. Real-time data assimilation in thermoacoustics was introduced by Nóvoa & Magri [10]. Statistical data assimilation methods are derived assuming that the model error is unbiased [11]. We generalise these techniques to low-order models, which typically provide a biased representation of the physics. Auto-regressive models have been employed to estimate the bias in modified Kalman Filters [12, 13], but they need to prescribe a-priori a model for the bias, and are limited to constant models or linear problems. We propose Echo State Networks (ESNs) to estimate the model bias. ESNs are recurrent neural networks based on reservoir computing [14, 15], which have been proved to learn complex temporal correlations in timeseries, such as nonlinear correlations, ergodic properties, or statistics of extreme events in turbulence [16–19]. They are general auto- regressive models [20]. This framework opens up new opportunities for exploiting experimental data for real-time state and parameter estimation of nonlinear thermoacoustics as well as for the estimation of model errors.

2. METHODOLOGY

Sequential data assimilation combines the prediction of a numerical model (the forecast ) with observable data to provide in real time a corrected state (the analysis ), which is a more accurate estimator of the physical state (the truth ). If the forecast is biased, e.g. if it is produced by a low-order model, the analysis remains in the biased subspace of the model. Here, we develop a method that eliminates the model bias. The biased forecast is mapped into an unbiased state, which is statistically combined with the observables to provide an unbiased analysis. This update indirectly determines the biased analysis, which is the initial condition of the next forecast step. Figure 1 conceptually illustrates this bias-aware sequential assimilation process.

Biased forecast

Biased analysis

Unbiased forecast Unbiased analysis

Truth

Observation d

State

State Time

Figure 1: Conceptual schematic of a bias-aware sequential data assimilation process. Truth (gray); observations (dark blue); biased forecast states (orange); biased analysis states (light blue); unbiased forecast states (red); and unbiased analysis states (purple). The circles represent pictorially the spread of the probability density functions: the larger the circles, the larger the uncertainty.

We propose a bias-aware data assimilation method to perform real-time inference of thermoacoustic state, parameters and model bias in a Rijke tube, which is a laboratory-scale horizontal open-ended tube with a heat source inside [10]. The observables are acoustic pressure measurements obtained with microphones [21, 22]. The ingredients for this framework are (i) a forecast low-order model, (ii) a statistical assimilation method, (iii) an estimator of the model bias, and (iv) data on acoustic pressure.

2.1. Forecast low-order model We use the dimensional form of the low-order model of a Rijke tube described in [23, 24]. The acoustic perturbations are assumed to evolve on top of a zero-mean flow, which is assumed to be a perfect gas with negligible e ff ect in the acoustics; the boundary conditions are assumed ideal; and we neglect viscous and body forces. The governing equations for the acoustic perturbations are the linearised momentum and energy conservation equations, which can be discretised in space using a Galerkin method [24]. We model the perturbation heat release rate as a compact source with the qualitative nonlinear time-delayed model

r u 0

3 − u ′ f ( t − τ )  − r u 0

˙ q ′ = β 

 , (1)

3

where β is the source strength, u 0 is the mean flow velocity, and u ′ f ( t − τ ) is the time delayed acoustic velocity at the flame location [25]. The forecast model can be written in compact notation as d ψ

d t = F ( ψ ) , ψ ( t = 0) = ψ 0 ,

y = M ( x ) ψ , (2)

where ψ ≡ ( η ; µ ; v ; α ) ∈ R N is the state vector, with N = 2 N m + N c + N α . This is a concatenation of the Galerkin acoustic modes for the velocity η ∈ R N m and pressure µ ∈ R N m ; v ∈ R N c , which are the modes resulting from the discretisation of an advection equation that numerically stores the delayed velocity [5]; and the thermoacoustic parameters α = ( β, τ ) ∈ R N α , which are the quantities that we aim to infer. The nonlinear operator is represented by F : R N → R N ; and y is the projection of the state vector into the observable space at the spatial locations x ∈ R N mic . The mapping is performed through the measurement operator matrix M ∈ R N mic × N defined though the Galerkin decomposition as

M = h 0 N mic × N m M N mic × N m 0 N mic × N c 0 N mic × N α i ,

with ( M ) ij = − sin ω j

! , i = 1 , ..., N mic , j = 1 , ..., N m , (3)

c 0 x i

where ω j is the frequency of each acoustic mode and c 0 the mean flow speed of sound.

2.2. Ensemble Square-Root Kalman Filter We employ the Ensemble Square-Root Kalman Filter (EnSRKF) as the statistical assimilation method, e.g. [10]. This is an ensemble-transform Kalman filter that is suitable for nonlinear problems, which avoids sampling errors of the observations [11, 26]. We take an ensemble approach to represent the state of the system, which reduces the computational cost compared to non-ensemble methods. For an ensemble A = ψ 1 , ψ 2 , . . . , ψ m
∈ R N × m , the expected state and its covariance are given by the sample statistics ψ ≈ 1 m P m j = 1 ψ j , and C ψψ ≈ 1 1 − m ΨΨ T , respectively, where Ψ =  ψ 1 − ψ , ψ 2 − ψ , . . . , ψ m − ψ  . The filter provides an analysis ensemble, A a , by updating

independently its mean, ψ a , and the deviations from the mean, Ψ a , as

ψ a = ψ f + Ψ f  M Ψ f  T  ( m − 1) C ϵϵ + M Ψ f  M Ψ f  T  − 1  d − y f  , (4a)

Ψ a = Ψ f V ( I − Σ ) 1 / 2 V T , (4b)

where the subscripts ‘f’ and ‘a’ stand for ‘forecast’ and ‘analysis’, respectively; d ∈ R N mic is the vector of the pressure observations; C ϵϵ ∈ R N mic × N mic is the observations’ error covariance matrix; y f is the unbiased expectation of the observables from the forecast model; the identity matrix is represented by I ; and V and Σ are the matrix of eigenvectors and of eigenvalues, from the singular value decomposition

V Σ V T =  M Ψ f  T  ( m − 1) C ϵϵ + M Ψ f  M Ψ f  T  − 1  M Ψ f  . (4c)

If the EnSRKF provides an analysis ensemble with non-physical parameters (e.g. negative time delays), we disregard the update and inflate the ensemble with a multiplicative inflation factor [10].

2.3. Echo State Network We define a time-dependent model bias, b ( t ) ∈ R N mic , as the di ff erence between the observations at each microphone, and the expected estimation from the low-order model. The echo state network predicts the evolution of this di ff erence by learning its temporal correlations. This is achieved by mapping at a time t the bias at the previous time step, b in ( t − ∆ t ), into a high-dimensional reservoir by the input matrix, W in ∈ R N r × ( N mic + 1) , where N r is the number of neurons in the reservoir. The state of the reservoir is defined by the matrix W ∈ R N r × N r , and a time-dependent high-dimensional vector r ( t ) ∈ R N r . The updated state vector is mapped back into the physical state by the output matrix, W out ∈ R ( N mic + 1) × N r . At a time t , the reservoir state and the bias are determined recursively as

r ( t ) = tanh  W in h ˜b ( t − ∆ t ); 0 . 1 i + Wr ( t − ∆ t )  and b ( t ) = W out [ r ( t ); 1] , (5)

where the operation [ ; ] represents vertical concatenation; ˜ b is the bias normalised by its range component–wise; and the constant 0.1 is used to break the symmetry of the network [27]. Wit h this definition of the bias, the expected unbiased forecast observables in Equation 4a are y f = M ψ f + b . The time-propagation of the bias can be performed in either open or closed loop (Figure 2). In open- loop, the input in each time step consists of the observed bias, while in closed-loop the estimated bias is used as the input for the next forecast step. During data assimilation, we only use the open-loop to get the ESN running, i.e., to have the state of the reservoir, r ( t ), up to date with the system (this procedure is known as washout). Then, sequential closed-loops are run in parallel with the EnSRKF, and are re-initialised with the measurement data.

(a) (b)

…

…

Figure 2: Schematic of the echo state network configurations: (a) open-loop and (b) closed-loop. The diamonds indicate inputs from the observable data.

The matrices W , W in , and W out are fixed during data assimilation. The state matrix, W , is a sparse and randomly generated Erd˝os-Renyi matrix. Each neuron, i.e., each row of W , has on average d = 5

Table 1: Summary of the parameters. Multiple values indicate that the parameters are optimised in the given range.

Data assimilation parameters ESN parameters

Ensemble size ( m ) 50 Training interval 1.2 s

State vector size ( N = 2 N m + N c + N α ) 20 + 50 + 2 Open-loop interval (washout) 0.025 s

Number of microphones ( N mic ) 6 Training β 5 · 10 6 Ws 1 / 2

m 5 / 2 Initial parameter uncertainty 15% Training τ 1 . 5 · 10 − 3 s

Observations standard deviation 1% Spectral radius ( ρ ) [0.5, 1.2]

Time between analysis 0.004 s Input scaling ( σ in ) [0.001, 2]

Inflation factor 1.01 Tikhonov param. (-log 10 γ t ) (4, 8, 12, 16)

non-zero elements sampled from a uniform distribution in [ − 1 , 1]. The state matrix is re-scaled by its spectral radius, ρ . Similarly, W in is sparse and randomly generated, but with only one non-zero element per row, sampled from a uniform distribution in [ − σ in , σ in ], where σ in is the input scaling. The output matrix W out is obtained from training, which consists of solving a linear system

( RR T + γ t I ) W T out = RB T out , (6)

for a training dataset of length N tr , where R ∈ R ( N r + 1) × N tr is the horizontal time-concatenation of the augmented reservoir state, i.e., [ r ( t i ); 1] with i = 1 , ..., N tr ; B out ∈ R N mic × N tr is the time-concatenation of the output data; and γ t is the Tikhonov regularisation parameter [14]. We add white noise with standard deviation, σ n = 0 . 03 σ u to the training data, where σ u is the component-wise standard deviation of the training data. This regularises further the network and improves the performance [14]. Finally, we optimise the parameters σ in , ρ and γ t with Recycle Validation [19]. Table 1 summarises the parameters used for the data assimilation and the echo state network.

2.4. Observables from a higher-order model We create the observable data from a travelling-wave model [28]. A detailed description of the implementation of this model can be found in [29]. The travelling-wave approach relaxes assumptions from the low-order model: (i) the mean flow e ff ects are not neglected, and jump conditions across the flame are implemented; (ii) the non-ideal boundary conditions are defined through reflection coe ffi cients; and (iii) the compact heat source assumption is relaxed by employing a flame kinematic model that kinematically tracks the flame front area [28].

3. RESULTS

We train the ESN with a 1.2 s dataset. The training data is the di ff erence between the statistically stationary limit cycle from the travelling–wave model, and a realisation of the lower–order model with an initial guess of the parameters β train and τ train . The training process optimises the weights of the output matrix as well as the network’s key parameters. To select the spectral radius, the input scaling and the Tikhonov parameter, we minimise the average mean-squared error (MSE) of multiple closed- loop predictions in the training data, which are defined through the Recycling Validation strategy [19] (Figure 3). The input scaling and spectral radius are tuned in the range [0 . 001 , 2] × [0 . 5 , 1 . 2]. The minimisation is performed with a Bayesian optimisation that starts from a 4 × 4 grid, and evaluates the function in four additional search points. In each point, the MSE error is computed for the

3.15 0

Input scaling (log-scale)

2.55

log 10 MSE

1

1.95

2

1.35

3

0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 0.75

Spectral Radius

Figure 3: Parameter optimisation through Recycle Validation and Bayesian optimisation. Contour plot of the mean-square error (MSE). The MSE is evaluated in a 4 × 4 grid (triangles), and four search points (squares).

Tikhonov parameters 10 − 4 , 10 − 8 , 10 − 12 , and 10 − 16 . The optimal parameter combination is found to be ρ = 0 . 9667 , σ in = 0 . 001 , and γ t = 1 · 10 − 16 . Figure 4 shows the results from state, parameter and model bias estimation. The ESN is initialised with an open-loop of 0.025 s sampled at 10000 Hz before the data assimilation begins (panel (a)). The biased solution is the expectation from the low-order model, and the unbiased filtered solution is the resulting pressure state after adding the model bias estimated with the ESN. The bias-aware data assimilation recovers the true pressure state and learns its dynamics. The attractor of the unbiased filtered solution, which is obtained after the filter reaches convergence, evolves around the true attractor, with some noise (panel (b)); and the dominant frequency spectrum is recovered (panel (d)). In addition, the parameters converges after 0.25 s to roughly the same time delay as the training value, and half of the training heat source strength (panel (c)). Therefore, the ESN is capable of estimating the bias of data that was not seen by the machine. (This unseen data has di ff erent thermoacoustic parameters and dynamics than the training dataset.) In the machine-learning terminology, the ESN generalises satisfactorily to unseen scenarios. This is an ESN’s key capability, which can be exploited in more realistic scenarios, for example, when the data comes from experiments.

4. CONCLUSIONS

The dominant dynamics of thermoacoustic oscillations can be modelled with qualitatively accurate low-order models. We propose a bias-aware data assimilation method to make reduced order models quantitatively correct. This is achieved by seamlessly combining an echo state network with a Bayesian sequential data assimilation method based on the ensemble Kalman filter. The model error (bias) is estimated by the echo state network, while the state-and-parameter estimation is given by the data assimilation algorithm, which combines the model and bias prediction with observations to provide the most likely estimation. The echo state network is trained a priori, then it is initialised before running in parallel with the sequential data assimilation algorithm. We test this framework with observations created from a higher-fidelity model, which includes mean flow e ff ects, non-ideal boundary conditions and a kinematic model of the flame. With a short training time, the echo state network is able to learn the dynamics of the physical information missing in the low-order model. The proposed bias-aware methodology successfully predicts, in real time, the time-accurate thermoacoustic state and its statistics, as well as the low-order model error and key thermoacoustic parameters. This work opens up possibilities for inferring model errors on the fly.

(b)

(a)

Truth Biased filtered sol. Unbiased filtered sol. Analysis step

Biased filtered sol. Truth Unbiased filtered sol.

5

6

p' ( t + 2 ζ ) . 10 -3

20

p' max ( i + 1) . 10 -3

5

p' mic 1 [Pa] . 10 -3

0

4

0

10

-5

1.4 1.42 1.44

2

-5

0

-5 0 5

-5 0 5

0 2 4 6

1 1.1 1.2 1.3 1.4 1.5 1.6 t [s]

p' max ( i ) . 10 -3 p' ( t ) . 10 -3 p' ( t + ζ ) . 10 -3

(c)

(d)

1.5

6000

β ( t ) / β train

t [1.3, 1.7] s

Truth Biased filtered sol. Unbiased filtered sol.

τ ( t ) / τ train

1.0

4000

PSD

2000

0.5

0

0

1 1.1 1.2 1.3 1.4 1.5 1.6 t [s]

Frequency [Hz] 0 50 100 150 200 250

Figure 4: State, parameter and model bias estimation. Time evolution of (a) the true (light grey), biased (dashed black) and unbiased (dark grey) filtered acoustic pressures at x / L = 0 . 18, the red circles are the observations, and the vertical dashed lines show the data assimilation window; and (c) the thermoacoustic parameters with their standard deviations, normalised to their training value. The (b) phase portrait and first return map, and (d) power spectral density (PSD) show the dynamics of the truth (dark blue), the unbiased (light blue) and biased (orange) filtered solutions after convergence.

REFERENCES

[1] Tim C. Lieuwen. Unsteady Combustor Physics . Cambridge University Press, Cambridge, 2012. [2] Lord Rayleigh. The Explanation of Certain Acoustical Phenomena. Nature , 18:319–321, July 1878. [3] Fred E. C. Culick. Unsteady motions in combustion chambers for propulsion systems. Technical report, NATO Research and Technology Organization Neuilly-Sur-Seine (France), 2006. [4] Lipika Kabiraj and R. I. Sujith. Nonlinear self-excited thermoacoustic oscillations: intermittency and flame blowout. Journal of Fluid Mechanics , 713:376–397, 2012. [5] Francisco Huhn and Luca Magri. Stability, sensitivity and optimisation of chaotic acoustic oscillations. Journal of Fluid Mechanics , 882:A24, 2020. [6] Luca Magri. Adjoint methods as design tools in thermoacoustics. Applied Mechanics Reviews , 71(2), March 2019. [7] Tullio Traverso and Luca Magri. Data Assimilation in a Nonlinear Time-Delayed Dynamical System with Lagrangian Optimization. Computational Science – ICCS 2019 , page 156–168, 2019. [8] Ross N. Bannister. A review of operational methods of variational and ensemble-variational data assimilation. Quarterly Journal of the Royal Meteorological Society , 143(703):607–633, 2017. [9] Jacqueline O’Connor, Vishal Acharya, and Timothy Lieuwen. Transverse combustion instabilities: Acoustic, fluid mechanic, and flame processes. Progress in Energy and Combustion

Science , 49:1–39, 2015. [10] Andrea Nóvoa and Luca Magri. Real-time thermoacoustic data assimilation. arXiv preprint arXiv:2106.06409 , 2021. [11] Geir Evensen. Data assimilation: the ensemble Kalman filter . Springer Science & Business Media, 2009. [12] Bernard Friedland. Treatment of bias in recursive filtering. IEEE Transactions on Automatic Control , 14(4):359–367, 1969. [13] Jean-Philippe Drécourt, Henrik Madsen, and Dan Rosbjerg. Bias aware kalman filters: Comparison and improvements. Advances in water resources , 29(5):707–718, 2006. [14] Mantas Lukoševiˇcius. A practical guide to applying echo state networks. In Neural networks: Tricks of the trade , pages 659–686. Springer, 2012. [15] Lyudmila Grigoryeva and Juan-Pablo Ortega. Echo state networks are universal. Neural Networks , 108:495–508, 2018. [16] Herbert Jaeger and Harald Haas. Harnessing nonlinearity: Predicting chaotic systems and saving energy in wireless communication. Science , 304(5667):78–80, 2004. [17] Wolfgang Maass, Thomas Natschläger, and Henry Markram. Real-time computing without stable states: A new framework for neural computation based on perturbations. Neural computation , 14(11):2531–2560, 2002. [18] Francisco Huhn and Luca Magri. Gradient-free optimization of chaotic acoustics with reservoir computing. Physical Review Fluids , 7(1):014402, 2022. [19] Alberto Racca and Luca Magri. Robust optimization and validation of echo state networks for learning chaotic dynamics. Neural Networks , 142:252–268, 2021. [20] Charu C Aggarwal et al. Neural networks and deep learning. Springer , 10:978–3, 2018. [21] Tim C. Lieuwen and Vigor Yang. Combustion instabilities in gas turbine engines: operational experience, fundamental mechanisms, and modelling . American Institute of Aeronautics and Astronautics, 2005. [22] Lipika Kabiraj, Aditya Saurabh, Pankaj Wahi, and R. I. Sujith. Route to chaos for combustion instability in ducted laminar premixed flames. Chaos: An Interdisciplinary Journal of Nonlinear Science , 22(2):023129, 2012. [23] Koushik Balasubramanian and R. I. Sujith. Thermoacoustic instability in a Rijke tube: Non- normality and nonlinearity. Physics of Fluids , 20(044103), April 2008. [24] Matthew P. Juniper. Triggering in the horizontal Rijke tube: non-normality, transient growth and bypass transition. Journal of Fluid Mechanics , 667:272–308, January 2011. [25] Maria A. Heckl. Non-linear Acoustic E ff ects in the Rijke Tube. Acustica , 72, 1990. [26] David M. Livings, Sarah L. Dance, and Nancy K. Nichols. Unbiased ensemble square root filters. Physica D: Nonlinear Phenomena , 237(8):1021–1028, 2008. [27] Francisco Huhn and Luca Magri. Learning ergodic averages in chaotic systems. In International Conference on Computational Science , pages 124–132. Springer, 2020. [28] Ann P. Dowling. A kinematic model of a ducted flame. Journal of fluid mechanics , 394:51–72, 1999. [29] José Guillermo Aguilar Pérez. Sensitivity analysis and optimization in low order thermoacoustic models . PhD thesis, University of Cambridge, 2019.