Sage Journals: Discover world-class research

Abstract

We present a method that combines multiple sensory modalities in a rocket thrust chamber to predict impending thermoacoustic instabilities with uncertainties. This is accomplished by training an autoregressive Bayesian neural network model that forecasts the future amplitude of the dynamic pressure time series, using multiple sensor measurements (injector pressure/ temperature measurements, static chamber pressure, high-frequency dynamic pressure measurements, high-frequency OH* chemiluminescence measurements) and future flow rate control signals as input. The method is validated using experimental data from a representative cryogenic research thrust chamber. The Bayesian nature of our algorithms allows us to work with a dataset whose size is restricted by the expense of each experimental run, without making overconfident extrapolations. We find that the networks are able to accurately forecast the evolution of the pressure amplitude and anticipate instability events on unseen experimental runs 500 milliseconds in advance. We compare the predictive accuracy of multiple models using different combinations of sensor inputs. We find that the high-frequency dynamic pressure signal is particularly informative. We also use the technique of integrated gradients to interpret the influence of different sensor inputs on the model prediction. The negative log-likelihood of data points in the test dataset indicates that prediction uncertainties are well-characterized by our model and simulating a sensor failure event results in a dramatic increase in the epistemic component of the uncertainty, as would be expected when a Bayesian method encounters unfamiliar, out-of-distribution inputs.

Keywords

Thermoacoustic instabilities multimodal data bayesian machine learning interpretable machine learning

Introduction

High-frequency thermoacoustic instabilities are a potentially catastrophic phenomenon in rocket and aircraft engines. They usually arise from the coupling between unsteady heat release rate and acoustic eigenmodes of the combustion chamber. Acoustic waves incident on a flame cause fluctuations in its heat release rate. If a flame releases more (less) heat than average during instants of higher (lower) local pressure, then more work is done by the gas during the acoustic expansion phase than is done on it during the acoustic compression phase. If this work is not dissipated then the oscillation amplitude grows and the system becomes thermoacoustically unstable.

Thermoacoustic instabilities can be challenging to model and predict because they are usually sensitive to small changes in a combustion system.¹ The large pressure fluctuations and accelerated heat transfer caused by these oscillations can lead to structural failure, because rocket engines have high energy densities and low factors of safety in their structural design. From the F1 rocket engine, which underwent costly redesigns in the 1960s² to the modern Japanese LE-9 engine,³ they continue to plague the development of liquid propellant rocket engines.

Designers can suppress thermoacoustic oscillations in two ways. The first method involves fitting passive devices such as acoustic liners,² baffles⁴ and resonators⁵ or the use of tuned passive control, where sensor measurements are used to diagnose thermoacoustic oscillations which can then be mitigated by adjusting a passive damping device or avoiding the operating point.⁶ The second method is active feedback control.⁷ The challenging demands of high power output, reliability and robustness placed on an active feedback control system means that feedback control is not used in rocket engines and aircraft engines. This paper proposes the diagnosis of impending instabilities using multiple sensor data-streams as input and Bayesian neural networks for forecasting. The forecast can be used by a rocket engine controller to perform a mass flow rate adjustment and avoid the unstable regions.

The remainder of the paper is paper is structured as follows. We provide an overview of related work in this field and explain the motivation behind the current study. Next, we describe the experimental rocket thrust chamber which generated our dataset. This is followed by details of the statistical tools we used: Bayesian neural network ensembles for forecasting and integrated gradients for interpreting them. We then compare the quality of the forecasts obtained using different model inputs and identify the influence of important input features on predictions using integrated gradients. We conclude with a brief summary of our findings and a discussion of potential avenues for future work.

The construction of instability precursors from pressure measurements and optical measurements has attracted considerable attention from the combustion research community, as reviewed by Juniper and Sujith.⁸ Lieuwen⁹ used the autocorrelation decay of combustion noise, filtered around an acoustic eigenfrequency, to obtain an effective damping coefficient for the corresponding instability mode. Other researchers have used tools from nonlinear time series analysis, which capture the transition from the chaotic behaviour displayed by stable turbulent combustors to the deterministic acoustics during instability, such as Gottwald’s 0-1 test¹⁰ and the Wayland test for nonlinear determinism¹¹ to assess the stability margin of a combustor. Nair et al.¹² reported that instability is often presaged by intermittent bursts of high-amplitude oscillations and used recurrence quantification analysis (RQA) to detect these. In a later paper, Nair and Sujith¹³ noted that combustion noise tends to lose its multifractality as the system transitions to instability. They proposed the Hurst exponent as an indicator of impending instability. Measures derived from symbolic time series analysis (STSA)¹⁴ and complex networks¹⁵ are also able to capture the onset of instability.

Recent studies have explored machine learning techniques to learn precursors of instability from data. These promise greater accuracy than physics-based precursors of instability, though at the cost of robustness and generalizability. Hidden Markov models constructed from the output of STSA¹⁶ or directly from pressure measurements¹⁷ have been used to classify the state of combustors. Hachijo et al.¹⁸ have projected pressure time series onto the entropy-complexity plane and used support vector machines (SVMs) to predict thermoacoustic instability. SVMs were also employed by Kobayashi et al.,¹⁹ who used them in combination with principal component analysis and ordinal pattern transition networks to build precursors from simultaneous pressure and chemiluminiscence measurements. Sengupta et al.²⁰ showed that the power spectrum of the noise can be used to predict the linear stability of a thermoacoustic eigenmode using Bayesian neural networks. Related work by McCartney et al.²¹ used the detrended fluctuation analysis (DFA) spectrum of the pressure signal as input to a random forest and found that this approach compared favorably to precursors from the literature. A recent study by Gangopadhyay et al.²² trained a 3D Convolutional Selective Autoencoder on flame videos to detect instabilities. Nóvoa and Magri²³ adopt a hybrid approach, using Bayesian data assimilation to learn the parameters of a low fidelity model and an echo state network to model the forecast bias. Transfer learning is also being explored by Mondal et al.²⁴ as a means of efficiently transferring knowledge of precursors across machines.

Most papers in the literature identify indicators of approaching instability. In this study, on the other hand, we use nonlinear autoregressive time series modeling with a Bayesian Neural Network to forecast, with uncertainties, the future amplitude of pressure fluctuations given the history of sensor signals and future flow rate control signals. This informs the engine controller of the expected timing and potential severity of a forthcoming instability. This approach also lets us take advantage of the rich instrumentation of our rig by integrating multimodal sensor data into our model: high frequency dynamic pressure ( $p^{'} (t)$ ) and OH* chemiluminescence ( $I^{'} (t)$ ) measurements, injector pressures ( $p_{H 2}$ , $p_{O 2}$ ) and temperatures ( $T_{H 2}$ , $T_{O 2}$ ), chamber static pressure ( $p_{cc}$ ) and fuel flow rates ( ${\dot{m}}_{O 2}, {\dot{m}}_{H 2}$ ). We compare models with different combinations of sensors, different sensor-derived features, and different lengths of signal histories in the inputs. We find that, although operating parameters and mass flow rate control signals can anticipate instabilities on their own, inclusion of dynamic pressure signal history, in particular, results in a measurable increase in the forecast accuracy. It is not possible to predict triggered thermoacoustic instabilities from operating parameters and control signals alone, so we need to integrate multimodal sensor data with operating parameters and control signals in an instability forecasting framework.

The size of our dataset is limited by the high cost and effort required for each additional test run on the rocket thrust chamber. Fortunately, the uncertainty-aware nature of the Bayesian Neural Network allows us to work in this medium-data regime without overconfident extrapolation. We compute log-likelihoods on the test dataset and find that the uncertainty is well-quantified. We also simulate a sensor failure to assess the robustness of our Bayesian forecasting tool. As expected, the network exhibits large epistemic uncertainties during the failure, signalling the unfamiliar out-of-distribution nature of the inputs. Finally we use the interpretation technique of integrated gradients (IG) to understand how the different input features in our model influence the predictions.

Experimental setup

The dataset used in this study is derived from experiments performed on the research combustor BKD²⁵ operated at the P8 test facility^26,27 of the DLR Institute of Space Propulsion in Lampoldshausen. It has three main components: an injector head, a cylindrical combustion chamber and a convergent-divergent nozzle, as shown in Figure 1. The combustion chamber is water-cooled and is designed to deal with the high thermal loads that are expected during instability events. The L42 injector head has 42 shear coaxial injectors and is operated with a liquid oxygen (LOX)/ hydrogen (LH2) propellant combination. The cylindrical combustion chamber is 80 mm in diameter and the nozzle throat diameter is 50 mm, resulting in a contraction ratio of 2.56. The chamber static pressure ( $p_{cc}$ ) is varied between $50 - 80$ bar and mixture ratio of oxidizer to fuel (ROF = ${\dot{m}}_{O 2} / {\dot{m}}_{H 2}$ ) between $2 - 6$ . The experiments considered here have a LOX injection temperature $T_{O 2}$ around 110 K and a hydrogen injection temperature $T_{H 2}$ around 100 K.

Figure 1.

Experimental thrust chamber BKD.

For the operating point with $p_{cc} = 80$ bar, ROF $= 6$ , the total propellant mass flow rate is 6.7 kg/s, the theoretical thermal power is 90 MW and the thrust achieved is about 24 kN. These specifications place BKD at the lower end of small upper stage engines.

A representative BKD test sequence, along with a spectrogram of dynamic pressure oscillations inside the combustion chamber, is shown in Figure 2. Stable and unstable operating conditions can be identified in the spectrogram. Strong high-frequency combustion instabilities of the first tangential (1T) mode at about 10 kHz were excited consistently when approaching the operating point at $p_{cc} =$ 80 bar, ROF $= 6$ . The coupling mechanism for this instability is injection-driven.²⁵ The flame dynamics are modulated by the LOX post acoustics and combustion instabilities emerge when the frequency of the 1T chamber mode matches the second longitudinal eigenmodes of the LOX posts. This mechanism has been confirmed using high-speed flame imaging.²⁸

Figure 2.

BKD test sequence and $p^{'}$ spectrogram showing self-excited instability of the first tangential mode.²⁵ $p_{cc}$ is the chamber static pressure, ROF is the mixture ratio of oxidizer to fuel, and $T_{H 2}, T_{O 2}$ are the LOX and hydrogen injection temperatures respectively.

The extreme conditions within the thrust chamber, where temperatures reach 3600 K and the pressure reaches 80 bar, limit the diagnostic instruments for our instability investigations. A specially designed measurement ring is placed between the injector head and the cylindrical combustion chamber segment, as shown in Figure 1. At this location, the temperatures are moderated by the injection of cryogenic propellants, which allows the mounted sensors to survive several test runs. Eight Kistler type 6043A water-cooled high-frequency piezoelectric pressure sensors are flush-mounted in the ring with an even circumferential distribution, in order to measure the chamber pressure oscillations $p^{'} (t)$ . The high-frequency pressure sensors have a measurement range set to $\pm 30$ bar and a sampling rate of 100 kHz. An anti-aliasing filter with a cutoff frequency of 30 kHz is applied. Three fibre-optical probes are used to record the OH* radiation intensity $I^{'} (t)$ of selected individual flames. The full acceptance angle of the optical probes is approximately $2^{\circ}$ . In order to capture the OH* radiation signal, the probes are equipped with inteference filters with a center wavelength of 310 nm. The sampling frequency of the $I^{'} (t)$ signals is also 100 kHz. There are also several low sampling frequency (100 Hz) sensors that measure the chamber static pressure $p_{cc}$ , the mass flow rates of fuel ${\dot{m}}_{H 2}$ and oxidizer ${\dot{m}}_{O 2}$ and injector temperatures $T_{H 2}, T_{O 2}$ .

Methods

Nonlinear autoregressive models with exogenous variables (NARX) for timeseries modeling

To forecast instabilities, we train NARX-style models²⁹ to capture the functional relationship between the history of the sensor data, control signals and future pressure fluctuations. These model the evolution of pressure fluctuations as follows:

log (p_{rms}^{'} (t_{m} + Δ t)) = f (x_{m}, x_{m - 1}, \dots x_{m - h}, u_{m + f}) + σ_{ϵ}

(1)Here,

p_{rms}^{'} (t_{m} + Δ t)

is the root mean square amplitude of the dynamic pressure fluctuation signal in a 50 ms window centered around

t_{m} + Δ t

. We choose

Δ t = 500

ms for our models because we consider this to be sufficient warning for our rocket engine controller to detect and react to an impending instability with a mass flow rate adjustment. The state variables

x_{m}, x_{m - 1}, \dots x_{m - h}

consist of various sensor data (or transformations thereof) at times

t_{m} - 50, t_{m} - 100, \dots t_{m} - 50 h

ms. The order of dependence

h

determines how far in the past we look. In this study, we explore models with

h = 0

and

h = 2

. The elements of

u_{m + f}

are the two mass flux control signals

{\dot{m}}_{H 2, s e t}

and

{\dot{m}}_{O 2, s e t}

averaged over the window between

t_{m}

and

t_{m} + Δ t

. These are the exogenous variables whose future values are known and can therefore be included in the inputs. It is possible to train models that capture higher order dependence on instantaneous values of the control signals instead of reducing them to an average, but we find that this does not improve predictive accuracy in our problem.

σ_{ϵ}

is the homoskedastic (constant) observation noise or aleatoric uncertainty term, which is assumed to be Gaussian. We forecast the logarithm of

p_{rms}^{'} (t_{m} + Δ t)

instead of

p_{rms}^{'} (t_{m} + Δ t)

itself because the noise in the pressure fluctuations increases with the amplitude, and so the simplifying assumption of constant homoskedastic noise only holds approximately once this logarithmic transform has been applied. It also forces

p_{rms}^{'}

predictions to be positive.

We model the nonlinear function $f$ in Equation (2) as a feedforward neural network with a simple multi-layer perceptron (MLP) architecture. In a standard neural network, one obtains a point estimate of the network parameters. Here, however, we train a machine learning algorithm on a dataset whose size is limited by the cost of performing experiments on large rocket thrust chambers. With a standard neural network this size limit could lead to naive and inaccurate extrapolations when the network encounters out-of-distribution inputs that are different from those which were experienced during training. In safety-critical systems such as rockets, control decisions based on overconfident predictions could be catastrophic. Bayesian Neural Networks (BNNs) are an elegant solution to this problem because they provide a confidence interval as well as a prediction.

Bayesian neural network ensembles

The Bayesian framework for training neural networks³⁰ involves placing a sensible prior probability distribution $p (θ)$ over the parameters $θ$ of the network. Observations $D$ are then used to update our belief about the parameter distribution. The posterior distribution $p (θ | D)$ is obtained using Bayes’ rule.

p (θ | D) = \frac{p (D | θ) p (θ)}{\int p (D | θ) p (θ) d θ}

(2)Bayesian models are resistant to overfitting on small datasets and are uncertainty aware. They are also useful for continual learning throughout the lifetime of a device in operation.³¹ Exact Bayesian inference, however, is not possible in neural network models because the integral

\int p (D | θ) p (θ) d θ

is analytically intractable. The most widely used numerical method to integrate over the posterior, Markov Chain Monte Carlo, is too computationally expensive to be practical. Variational inference is often used,³² in which the posterior is approximated using a convenient parametrization and the Kullback-Leibler divergence between this variational distribution and the true posterior is minimized using backpropagation. However, although computationally cheap, mean-field variational inference also has drawbacks such as not capturing correlations between parameters. Recently, a different approximate inference method, based on ensembling, has been proposed. This is called the anchored ensembling algorithm. It is cheap, simple and scalable but manages to outperform variational inference in several uncertainty quantification benchmarks.³³ Consider a dataset of

N_{D}

data points

(x_{n}, y_{n})

, where each data point consists of features

x_{n} \in R^{D}

and output

y_{n} \in R

. Define the likelihood for each data point as

p (y_{n} ∣ θ, x_{n}, σ_{ϵ}^{2}) = N (y_{n} ∣ NN (x_{n}; θ), σ_{ϵ}^{2})

, where

NN

is a neural network whose weights and biases form the latent variables

θ

while

σ_{ϵ}^{2}

is the observation noise. Define the prior on the weights and biases

θ

to be the standard normal

p (θ) = N (θ ∣ μ_{p r i o r}, Σ_{p r i o r}) .

The anchored ensembling algorithm (Figure 3) then does the following:

The parameters $θ_{a n c, j}$ of each $j$ -th member of our neural network ensemble are initialized by drawing from the prior distribution $N (μ_{p r i o r}, Σ_{p r i o r})$ .

Each ensemble member is trained using a modified loss function that anchors the parameters to their initial values. The loss function for the $j$ -th ensemble member is given by $L o s s_{j} = \sum_{i = 1}^{N_{D}} (y_{i} - {\hat{y}}_{j} (x_{i}, t_{i}))^{2} + ‖ Σ_{p r i o r}^{- 1 / 2} (θ_{j} - θ_{a n c, j}) ‖_{2}^{2}$ , where ${\hat{y}}_{j} (x_{i}, t_{i})$ is the neural network output and the $i$ -th diagonal element of $Σ$ is the ratio of data noise to the prior variance of the $i$ -th parameter.

Figure 3.

Approximate Bayesian inference using ensembles of neural networks.

Pearce et al.³³ prove that this procedure approximates the true posterior distribution for wide neural networks and that the trained neural networks in the ensemble may be treated as samples from an approximate posterior distribution. The resulting ensemble has predictions that converge when they are well-supported by the training data and diverge when they are not. The variance of the ensemble predictions is an estimate of the Bayesian model’s epistemic uncertainty. Epistemic or knowledge uncertainty is the reducible portion of uncertainty which stems from uncertainty in model parameters and decreases as more data is added. The total uncertainty of each prediction is obtained by adding the epistemic uncertainty and the irreducible aleatoric uncertainty $σ_{ϵ}$ which stems from observation noise:

σ_{total}^{2} = σ_{ϵ}^{2} + \frac{1}{N_{e}} \sum_{j} {\hat{y}}_{j}^{2} - (\frac{1}{N_{e}} \sum_{j} {\hat{y}}_{j})^{2}

(3)

Interpretation using integrated gradients

Neural networks have a reputation as black boxes, which disincentivizes their application to cases in which practitioners must understand why the algorithm is or is not working. We use the technique of integrated gradients (IG)³⁴ to attribute the predictions of our network ensemble to the input features. This is a simple scalable method that only requires access to the gradient operation. IG has several desirable properties that other attribution methods lack, such as implementation invariance, sensitivity, linearity, completeness and preservation of symmetry.

For deep neural networks, the gradients of the output with respect to the input are an analog of the linear regression coefficients. In a linear model, the regression coefficients characterize the relationship between input and output variables globally. In a nonlinear neural network, however, the gradient of the output at a point merely characterizes the local relationship between a predictor variable and the output. IG computes the path integral of the gradient of the outputs with respect to the inputs around a baseline to the input under consideration. For an image recognition algorithm, a completely black image could be a reasonable choice of baseline, while for a regression problem like ours, where the input variables were normalized to have zero mean and unit standard deviation, the average input of all zeros is a sensible baseline. We consider the straight-line path (in the feature space $R^{p}$ ) from the baseline $x^{'}$ to the input $x$ being considered. The IG attribution $a t t r_{j}$ for the $j$ -th feature is then defined as follows:

a t t r_{j} (x) = (x_{j} - x_{j}^{'}) \int_{α = 0}^{1} \frac{\partial f (x^{'} + α (x - x^{'}))}{\partial x_{j}} d α

(4)This integral is computed numerically by sampling the gradient evenly along the path. For our neural network ensemble, we average the integrated gradient attributions obtained from individual ensemble members to obtain a mean

a t t r_{j}

for each feature.

Transforming sensor signals into neural network inputs

We compare NARX models with different combinations and transformations of the sensor data in their state variables $x_{m - i}$ (Equation 2) and evaluate their accuracy.

Models that include the low frequency signals $p_{cc}, {\dot{m}}_{H 2}, {\dot{m}}_{O 2}, T_{H 2}$ and $T_{O 2}$ in their $x_{m - i}$ vectors use the average value of the signal between $t_{m} - 50 i$ and $t_{m} - 50 (i + 1)$ . For the high frequency pressure und OH* signals, $p^{'}$ and $I^{'}$ , we use transformations from the literature that are known to capture instability-relevant information as described next. The current $p_{rms}^{'}$ is a natural candidate which is included in all models that incorporate pressure information. Sengupta et al.²⁰ and McCartney et al.³⁵ have shown that, before a combustor transitions to instability, the power spectral density shows measurable changes which can be used as a precursor. We estimate the spectral density of the $p^{'}$ ( $p_{FFT}^{'}$ ) and $I^{'}$ ( $I_{FFT}^{'}$ ) time series by applying Welch’s method³⁶ to a segment of the signals between $t_{m} - 50 i$ ms and $t_{m} - 50 (i + 1)$ ms. In this study, a Hann window of length 30 with 50% overlap is used for both $p^{'}$ and $I^{'}$ .

Kobayashi et al.¹⁹ applied machine learning to the ordinal partition transition network of the $p^{'}$ and $I^{'}$ timeseries to detect instabilities early. The ordinal partition transition networks for an $m$ -dimensional time series are expressed by a weighted adjacency matrix $W$ consisting of $W_{i j} = P (π_{i} ⟶ π_{j})$ , $i, j \in [1, 2^{m}]$ , $\sum W_{i j} = 1$ , where $P (π_{i} ⟶ π_{j})$ is the observed probability from the $i$ th- to $j$ th-order transition patterns. In this study, the order patterns are $π_{1}$ , $π_{2}$ , $π_{3}$ , and $π_{4}$ which account for all combinations of the signs of the increments for our two-dimensional time series $[p^{'} (t), I^{'} (t)]$ . The probability distribution of the corresponding 16-dimensional transition patterns $W$ computed from the segment of the signals between $t_{m} - 50 i$ and $t_{m} - 50 (i + 1)$ are used in the state variable vector.

We also include a naive baseline model $p_{rms}^{'} (t_{m} + Δ t) = p_{rms}^{'} (t_{m})$ for comparison. This model, by design, is incapable of forecasting any future instabilities but, because large portions of the run consist of steady operation, its predictions have a reasonably good root mean squared error (RMSE). This model exists to provide a baseline RMSE and any useful forecast needs to be significantly more accurate than this.

Train-test-validation split, hyperparameters, performance metrics

We use a dataset of five experimental runs: three that are 50 seconds long and two that are 90 seconds long. Each run has two occurrences of the 1T instability, where an instability event is defined as being when the peak-to-peak amplitude of the pressure fluctuations exceeds 5% of the chamber static pressure.

To avoid data leakage, one experimental run is exclusively reserved for hyperparameter tuning. The timeseries is split into contiguous 250 millisecond blocks and every fifth block is used to evaluate the negative log-likelihood, which is the metric we use for tuning the BNN hyperparameters. The hyperparameters are the network depth, layer width, noise $σ_{ϵ}$ , the learning rate for the ADAM optimizer, the number of epochs and the number of neural networks in the ensemble. For simplicity, and to facilitate comparison between models, the same hyperparameters are chosen for all models. This does not compromise the performance of any particular model, since increasing the key hyperparameters, width and depth, benefits all models. We choose a 3 hidden layer MLP with 100 units in each layer because larger networks show only minimal gains in performance. A learning rate of $2 \times 10^{- 4}$ is found to be optimal for convergence. The aleatoric noise parameter $σ_{ϵ}$ is set to 0.07. In anchored ensembling, we train all members of the ensemble until convergence so the number of epochs was set to 512, which is found sufficient for all models. Converged estimates of the test data log-likelihood are obtained using 25 neural networks per ensemble.

To evaluate model performance, leave-one-out cross validation (LOOCV) is performed on the remaining four runs. Time series data involves strong temporal correlations, so assigning training and test data points randomly is dishonest because then the two sets would be highly similar. The test set must therefore be a completely independent experimental run. Additionally, LOOCV is necessary because with only four runs, firm conclusions cannot be drawn based on performance metrics for one particular test-train split, since any claimed efficacy or differences between models could be caused by chance. We report average performance metrics across the four possible splits.

To evaluate model accuracy, we compute the root mean squared error (RMSE) of our pressure amplitude forecast for the entire test run (full RMSE). We also calculate this on a subset of the test run comprising the two segments of the timeseries one second before and after a transition to instability event (transition RMSE). The RMSE is computed separately on this subset because our primary interest is in forecasting instabilities and the relative merits of different models become clearer on this subset where there is a dramatic change in the amplitude. As a measure of the quality of uncertainty, we report mean negative log-likelihood per test data point.

Results

The models were trained on a laptop with 16 GB RAM, an Intel i7-10870H CPU and an NVIDIA RTX 2070 GPU. Each NARX model took $\sim 2.5$ hours to train. The inference time of the ensemble averaged around $10$ milliseconds per sample, indicating that the algorithm is suitable as a real-time diagnostic tool, once it has been trained on experimental data.

Table 1 shows the average prediction RMSEs for the NARX models computed during cross-validation across the four runs. We note, firstly, that the NARX models outperform the naive baseline, which means that the data contains a predictable signal that can be extracted. We also observe that the simple models, which use only the two operating point parameters $p_{cc}$ and ROF as inputs, perform well. This is unsurprising because we know that the 1T instabilities are consistently triggered when the operating point at $p_{cc} = 80$ bar and ROF $= 6$ is approached. This means that knowing the current operating point and future flow rate control signals is informative. All transitions to instability in the test runs differ from each other, however, even though they have identical control trajectories. The cross-validation shows unambiguously that the data from the additional instrumentation improves the prediction of future pressure fluctuations. Models that include the dynamic pressure spectrum $p_{FFT}^{'}$ , radiation intensity spectrum $I_{FFT}^{'}$ , current amplitude $p_{rms}^{'}$ , ordinal network transition probabilities, $W$ , $T_{O 2}$ , $T_{H 2}$ , $p_{O 2}$ or $p_{H 2}$ , have better prediction accuracies, especially on the transition subset. Including higher order dependencies on state variables also provides a slight boost in accuracy. Interestingly, the dynamic pressure spectrum appears more informative than the spectrum of the OH* radiation intensity. This could be due to the fact that a single optical probe observes only a small volume of the chamber, while the pressure field integrates information across the whole chamber.

Table 1.

Cross-validation RMSEs (in bar) for entire runs (Full RMSE) and the transition subset of the run (Transition RMSE) and mean negative log-likelihoods (NLL) for different NARX models. Lower RMSEs and NLLs are better.

State variable inputs $x$	Order	Full RMSE	Transition RMSE	NLL
Baseline	–	0.1552	0.2134	–
$p_{cc}$ , ROF	0	0.1289	0.1630	$-$ 1.622
$p_{cc}$ , ROF	2	0.1202	0.1576	$-$ 1.631
$p_{rms}^{'}$ , $p_{FFT}^{'}$	0	0.1197	0.1546	$-$ 1.5752
$p_{rms}^{'}$ , $p_{FFT}^{'}$	2	0.1169	0.1508	$-$ 1.5467
$p_{rms}^{'}$ , $I_{FFT}^{'}$	0	0.1255	0.1612	$-$ 1.5154
$p_{rms}^{'}$ , $W$	0	0.1241	0.1566	$-$ 1.5801
$T_{O 2}, T_{H 2}, p_{O 2}, p_{H 2}$	0	0.1243	0.1578	$-$ 1.5948
$p_{cc}$ , ROF, $p_{rms}^{'}$ , $p_{FFT}^{'}$ , $T_{O 2}, T_{H 2}, p_{O 2}, p_{H 2}$	0	0.1048	0.1407	$-$ 1.5632
$p_{cc}$ , ROF, $p_{rms}^{'}$ , $p_{FFT}^{'}$ , $T_{O 2}, T_{H 2}, p_{O 2}, p_{H 2}$	2	0.1005	0.1311	$-$ 1.5350

Table 1 also reports the mean negative log-likelihood (NLL) per datapoint, averaged across the four cross-validations. The mean negative log-likelihoods are fairly low for most of our models, implying that very few observations lie far outside the uncertainty bounds of the forecast (Figure 4). This is also confirmed by the percentage of data points in the test set within the $\pm 1$ standard deviation (s.d.) and $\pm 3$ s.d. bounds, which are roughly in line with our Gaussian assumptions (68% within 1 s.d., 99.7% within 3 s.d.). Additionally, though models using more features have higher prediction accuracies, they have slightly worse NLLs in some cases because the presence of additional features results in higher epistemic uncertainties.

Figure 4.

Forecast $p_{rms}^{'}$ values for the four experimental runs based on sensor data prior to $t - 0.5$ sec, superimposed on observed $p_{rms}^{'}$ values at time $t$ . The model was trained on the other three runs in each case. This model used $p_{FFT}^{'}$ and $p_{rms}^{'}$ as state variables and order $h = 2$ . Actual observations are shown as black dots, 3 standard deviation (S.D.) total uncertainty bounds for the forecast are in light blue and 3 standard deviation (S.D.) epistemic uncertainty bounds are in deep blue.

To test the robustness of the Bayesian model against out-of-distribution inputs, we simulate a sensor failure event. This is shown in Figure 5. To simulate an optical probe icing event between 10 sec and 40 sec, the original photomultiplier signal is replaced by white noise with mean 0.05 and amplitude 0.015. During this simulated probe failure, the epistemic uncertainty of the Bayesian neural network greatly increases, indicating that the inputs are highly dissimilar to conditions experienced during training. Predictions remain reasonable despite the uninformative optical signal, although they become much less accurate.

Figure 5.

Bayesian Neural Network epistemic uncertainty (top) and $p_{rms}^{'}$ forecasts (bottom) for experimental run no. 4 where a fibre optical probe failure was simulated between $t = 10$ sec and $t = 40$ sec (marked by two grey dotted lines). The model used $I_{FFT}^{'}, p_{rms}^{'}$ as input state variables and order = 0. Actual observations are shown as black dots, 3 standard deviation (S.D.) total uncertainty bounds for the forecast are represented by light blue bars and 3 standard deviation (S.D.) epistemic uncertainty bounds are in deep blue. The light blue bars are barely discernible from the deep blue bars in the failure region because the total uncertainty is almost entirely dominated by the large epistemic uncertainty.

We use the technique of integrated gradients to produce feature-level attribution plots, which can tell us how much a particular predictor influenced the prediction for a particular input. For example, the attribution plots in Figure 6 shows this technique being applied to two data points. The first one is a data point just prior to the first instability event in a run, where a large jump in the dynamic pressure fluctuation amplitude is forecast by the model. This model uses $p_{FFT}^{'}$ , $p_{rms}^{'}$ and the flowrate control signals as input. The most prominent feature is the large positive attribution for the oxygen flowrate control signal ${\dot{m}}_{O 2}$ , whose feature value is also large for this datapoint. This is expected, because a large ${\dot{m}}_{O 2}$ would increase the $p_{cc}$ and ROF, bringing the combustor closer to instability. The prediction is also modulated by contributions from the $p_{FFT}^{'}$ features. This datapoint has a concentration of power around the 10 kHz spectrum and low power in the mid to higher frequencies. There is a positive contribution from the 10 kHz frequency component of the spectrum, which is the frequency of the 1T mode, as well as from the 16.66k Hz frequency. This shows that the model has learned that an increased concentration of acoustic power around the frequency of the 1T mode and away from the mid-frequencies is an indication that the system is close to instability. The second data point is from stable operation. The attributions in this case have mostly negative values because the prediction of the model is strongly negative. The current pressure fluctuation amplitude $p_{rms}^{'}$ is close to $- 1.0$ in this case and naturally has a large negative attribution. The ${\dot{m}}_{H 2}$ input (here $- 0.78$ ) also contributes negatively because low ROFs are stable. Finally, the high concentration of power in the mid-frequencies also tells the network that this operating point will stay stable.

Figure 6.

Integrated gradient attributions and input feature values (normalized to have zero mean and unit standard deviation) for two datapoints from experimental run no. 4 at t = 16.0 seconds (top, immediately preceding an instability) and t = 60.0 seconds (bottom, stable). The model used $p_{FFT}^{'}$ and $p_{rms}^{'}$ as state variables and order = 0. The blue dots show the normalized values of the features for this particular datapoint and the bars represent the integrated gradient attributions for each feature. For example, feature 6 (spectral power around the 16.66 kHz frequency) in the top plot has a large negative value and a positive integrated gradients attribution, meaning that the network has learned that low power in this frequency band is a potential indicator of instability.

Conclusions

In this paper, we present the application of an uncertainty-aware machine learning algorithm to forecast instabilities using multimodal sensor data in a realistic, thermoacoustically unstable rocket thrust chamber. We predict the amplitude of dynamic pressure fluctuations $500$ ms in advance using temperature, static pressure, high-frequency dynamic pressure and OH* radiation intensity data recorded during test runs performed with the cryogenic rocket thrust chamber BKD. We perform a rigorous cross-validation of the autoregressive Bayesian neural network and find that our models are able to anticipate all the instability events with varying accuracy ( $25 - 40 %$ lower RMSEs on the transition subset compared to baseline), depending on the combination of inputs chosen. The inclusion of longer histories and additional sensor information, in general, boosts predictive accuracy. Models that use the dynamic pressure signal, in particular, show a measurable improvement in the forecast, confirming the findings of the literature on this topic. We gain insight into the predictions made by the neural network by using Integrated Gradients, which reveals the contribution made by each input feature to the final prediction.

A key benefit of our Bayesian approach is the ability to quantify the uncertainty in our predictions. In safety-critical machines such as rocket engines, a diagnostic tool must be robust to unexpected events such as sensor failures. We compute the log-likelihood of the test data given model predictions and find that uncertainties are well-characterized. We also simulate a sensor failure event in a model whose inputs included the spectrum of the fibre-optical probe signal and demonstrate the robustness of Bayesian algorithms. Integrated Gradients helps us identify important features and how they influence a particular prediction of the neural network.

Future work will focus on exploring the effectiveness of these tools on DLR’s other experimental datasets, such as the LOX-methane tests, which also displayed injection-coupled acoustic instabilities.³⁷ We will also explore the generalizability of data-driven diagnostics across multiple datasets and combustors, using approaches that utilise transfer learning or domain invariance. Hardware upgrades are currently being undertaken at the test facility that will allow the implementation of these methods into the engine control algorithm.³⁸

Footnotes

Acknowledgements

The authors would like to thank the crew of the P8 test bench and also Stefan Gröning for preparing and performing the test runs. Furthermore, we thank Wolfgang Armbruster for helping us process and understand the data. Ushnish Sengupta is an Early Stage Researcher within the MAGISTER consortium which receives funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 766264.

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iDs

Ushnish Sengupta

Justin Hardi

References

Juniper

. Sensitivity analysis of thermoacoustic instability with adjoint helmholtz solvers. Phys Rev Fluids 2018; 3: 110509.

Oefelein

Yang

. Comprehensive review of liquid-propellant combustion instabilities in F-1 engines. J Propuls Power 1993; 9: 657–677.

Watanabe

Manako

Onga

, et al. Combustion stability improvement of le-9 engine for booster stage of h3 launch vehicle. Mitsubishi Heavy Ind Tech Rev 2016; 53: 28.

Eldredge

Dowling

. The absorption of axial acoustic waves by a perforated liner with bias flow. J Fluid Mech 2003; 485: 307.

Zinn

. A theoretical study of nonlinear damping by Helmholtz resonators. J Sound Vib 1970; 13: 347–356.

Kobayashi

Murayama

Hachijo

, et al. Early detection of thermoacoustic combustion instability using a methodology combining complex networks and machine learning. Phys Rev Appl 2019; 11: 064034.

Dowling

Morgans

. Feedback control of combustion oscillations. Annu Rev Fluid Mech 2005; 37: 151–182.

Juniper

Sujith

. Sensitivity and nonlinearity of thermoacoustic oscillations. Annu Rev Fluid Mech 2018. DOI: 10.1146/annurev-fluid-122316-045125.

Lieuwen

. Online combustor stability margin assessment using dynamic pressure data. J Eng Gas Turbine Power 2005SEP. DOI: 10.1115/1.1850493.

10.

Nair

Thampi

Karuppusamy

, et al. Loss of chaos in combustion noise as a precursor of impending combustion instability. Int J Spray Combus Dyn 2013; 5: 273–290.

11.

Gotoda

Nikimoto

Miyano

, et al. Dynamic properties of combustion instability in a lean premixed gas-turbine combustor. Chaos 2011. DOI: 10.1063/1.3563577.

12.

Nair

Thampi

Sujith

. Intermittency route to thermoacoustic instability in turbulent combustors. J Fluid Mech 2014; 756: 470.

13.

Nair

Sujith

. Multifractality in combustion noise: predicting an impending combustion instability. J Fluid Mech 2014; 747: 635–655.

14.

Sarkar

Chakravarthy

Ramanan

, et al. Dynamic data-driven prediction of instability in a swirl-stabilized combustor. Int J Spray Combus Dyn 2016; 8: 235–253.

15.

Murugesan

Sujith

. Detecting the onset of an impending thermoacoustic instability using complex networks. J Propuls Power 2016; 32: 707–712.

16.

Jha

Virani

Ray

. Markov modeling of time series via spectral analysis for detection of combustion instabilities. In: Handbook of Dynamic Data Driven Applications Systems. Springer, 2018, pp. 123–138.

17.

Mondal

Ghalyan

Ray

, et al. Early detection of thermoacoustic instabilities using hidden Markov models. Combust Sci Technol 2018.

18.

Hachijo

Masuda

Kurosaka

, et al. Early detection of thermoacoustic combustion oscillations using a methodology combining statistical complexity and machine learning. Chaos: Interdiscipl J Nonlinear Sci 2019; 29: 103123.

19.

Kobayashi

Murayama

Hachijo

, et al. Early detection of thermoacoustic combustion instability using a methodology combining complex networks and machine learning. Phys Rev Appl 2019; 11: 064034.

20.

Sengupta

Rasmussen

Juniper

. Bayesian machine learning for the prognosis of combustion instabilities from noise. J Eng Gas Turbine Power 2020; 143: 071001.

21.

McCartney

Indlekofer

Polifke

. Online prediction of combustion instabilities using machine learning. J Eng Gas Turbine Power 2021; 143(7): 071001.

22.

Gangopadhyay

Ramanan

K Boor

, et al. 3d convolutional selective autoencoder for instability detection in combustion systems. Energy and AI 2021; 4: 100067.

23.

Nóvoa

Magri

. Real-time thermoacoustic data assimilation. J Fluid Mech 2022; 948: A35.

24.

Mondal

Chattopadhyay

Mukhopadhyay

, et al. Transfer learning of deep neural networks for predicting thermoacoustic instabilities in combustion systems. Energy and AI 2021; 100085. DOI: 10.1016/j.egyai.2021.100085. https://www.sciencedirect.com/science/article/pii/S2666546821.

25.

Gröning

Hardi

Suslov

, et al. Injector-driven combustion instabilities in a hydrogen/Oxygen rocket combustor. J Propuls Power 2016; 32: 560–573.

26.

Fröhlke

Haberzettl

Haidn

et al. First hot fire test campaign at the french/german research facility p8. In: 33rd Joint Propulsion Conference and Exhibit. 1997, p. 2929.

27.

Koschel

Haidn

Krülle

. P8-the new French/German test facility for H2/O2-high pressure rocket engine combustion research. In: 11th World Hydrogen Energy Conference, Stuttgart, June 23-28, 1996. DECHEMA, Frankfurt, Germany, 1996, pp. 1815–1825.

28.

Armbruster

Hardi

Suslov

, et al. Injector-driven flame dynamics in a high-pressure multi-element oxygen-hydrogen rocket thrust chamber. J Propuls Power 2019; 35: 632–644.

29.

Lin

Horne

Tino

, et al. Learning long-term dependencies in narx recurrent neural networks. IEEE Trans Neural Netw 1996; 7: 1329–1338.

30.

Neal

. Bayesian Learning for Neural Networks, Lecture Notes in Statistics, volume 118. New York, NY: Springer New York, 1996. ISBN 978-0-387-94724-2. doi:10.1007/978-1-4612-0745-0. http://link.springer.com/10.1007/978-1-4612-0745-0.

31.

Barnaghi

Enshaeifar

, et al. Continual learning using Bayesian neural networks, 2019.

32.

Blundell

Cornebise

Kavukcuoglu

, et al. Weight Uncertainty in Neural Networks; 2015. http://arxiv.org/abs/1505.05424.

33.

Pearce

Zaki

Brintrup

, et al. Uncertainty in Neural Networks: Bayesian Ensembling, 2018. http://arxiv.org/abs/1810.05546.

34.

Sundararajan

Taly

Yan

. Axiomatic attribution for deep networks. In: International Conference on Machine Learning. PMLR, 2017, pp. 3319–3328.

35.

McCartney

Sengupta

Juniper

. Reducing uncertainty in the onset of combustion instabilities using dynamic pressure information and bayesian neural networks. In: Proceedings of ASME Turbo Expo 2021. ASME, 2021.

36.

Welch

. The use of fast fourier transform for the estimation of power spectra: a method based on time averaging over short, modified periodograms. IEEE Trans Audio Electroacoust 1967; 15: 70–73.

37.

Klein

Börner

Hardi

, et al. Injector-coupled thermoacoustic instabilities in an experimental LOx-methane rocket combustor during start-up. CEAS Space Journal 2020; 1–13.

38.

Waxenegger-Wilfing

Dresia

Oschwald

, et al. Hardware-in-the-loop tests of complex control software for rocket propulsion systems. In: 71st International Astronautical Congress, Virtual. IAC, 2020.

Forecasting thermoacoustic instabilities in liquid propellant rocket engines using multimodal Bayesian deep learning

Abstract

Keywords

Introduction

Experimental setup

Methods

Nonlinear autoregressive models with exogenous variables (NARX) for timeseries modeling

Bayesian neural network ensembles

Interpretation using integrated gradients

Transforming sensor signals into neural network inputs

Train-test-validation split, hyperparameters, performance metrics

Results

Conclusions

Footnotes

Acknowledgements

Declaration of Conflicting Interests

Funding

ORCID iDs

References