Early Detection of Thermoacoustic Instabilities in a Cryogenic Rocket Thrust Chamber using Combustion Noise Features and Machine Learning

11/25/2020 ∙ by Günther Waxenegger-Wilfing, et al. ∙ DLR 0

Combustion instabilities are particularly problematic for rocket thrust chambers because of their high energy release rates and their operation close to the structural limits. In the last decades, progress has been made in predicting high amplitude combustion instabilities but still, no reliable prediction ability is given. Reliable early warning signals are the main requirement for active combustion control systems. In this paper, we present a data-driven method for the early detection of thermoacoustic instabilities. Recurrence quantification analysis is used to calculate characteristic combustion features from short-length time series of dynamic pressure sensor data. Features like the recurrence rate are used to train support vector machines to detect the onset of an instability a few hundred milliseconds in advance. The performance of the proposed method is investigated on experimental data from a representative LOX/H_2 research thrust chamber. In most cases, the method is able to timely predict two types of thermoacoustic instabilities on test data not used for training. The results are compared with state-of-the-art early warning indicators.



There are no comments yet.


page 3

page 4

page 8

page 9

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

I Introduction

Future space missions and their increasing complexity require rocket engines with the capability of adjustment to the specific missions and therefore a wide envelope of operating conditions [6]. Today rocket engines are built and tested in terms of operational stability for the operating conditions they are explicitly built for. Due to the complexity of upcoming space missions reusing an engine without extensive testing for the new requirements which have to be met can result in a decisive advantage.

The occurrence of combustion instabilities presents high technical risks during rocket engine development projects and their operation [19]. Combustion instability refers to high-amplitude pressure oscillations during a combustion process and is an issue for all chemical propulsion systems, e.g. gas turbines. Combustion instabilities are particularly problematic for rocket thrust chambers due to their operation close to the structural limits caused by the necessity of extreme lightweight construction [20]. Those instabilities can be further sub-categorized into high-frequency (screeching) and low-frequency instabilities [19, 52]. In the latter case, a further distinction in chugging and pogo instabilities is applied. Pogo is characterized by a change of the external loads (e.g. thrust modulations) which are transferred to the structure and chugging is normally attributed to an interaction between the combustion process and propellant feed system elasticity. In the case of high-frequency instabilities, the frequency of pressure oscillations is typically linked to the acoustic modes of the combustion chamber. The existing feedback loop within these acoustic oscillations was first described by Lord Rayleigh [38]. If the pressure oscillation and the heat release are in phase the necessary condition for growing amplitudes is given. Otherwise, a damping effect applies. Due to the enormous power density in combustion devices for rockets, the transfer of a small amount of the total heat release rate into the acoustic field can cause rapidly growing pressure oscillations to damagingly high amplitudes. In the last decades, progress has been made in predicting high amplitude combustion instabilities but still, no absolute certainty has yet been achieved [41].

Two different approaches can be used to prevent the occurrence of thermoacoustic instability. On the one hand, the development of these oscillations can be suppressed via passive control (for example baffles or acoustic liners [19]) or on the other hand through active control. In the case of active control, the system has to be able to process the data in time by using available sensor time-series data. Questions regarding suitable sensors, actuators, and optimal control targets must also be clarified Furthermore, high robustness and reliability are mandatory. Due to the challenging requirements for an active control system no usage of such a system for an engine in operational service is known to the authors.

The construction of instability precursors from pressure measurements has been intensively studied by the thermoacoustic instability community [22]. The first step in this direction was taken by Lieuwen [26], who used the autocorrelation decay of combustion noise filtered around an acoustic eigenfreuqency to obtain an effective damping coefficient that indicated closeness to instability. Other researchers have borrowed measures from nonlinear time series analysis which capture the transition from the chaotic behaviour displayed by stable turbulent combustors to the deterministic acoustics during instability, e.g. Gottwald’s 0-1 test [33]. Similarly, Gotoda and coworkers have used the Wayland test for nonlinear determinism [9] to assess the stability margin of a combustor. Nair et al [34] reported that instability is often presaged by intermittent bursts of high-amplitude oscillations and used recurrence quantification analysis to detect these. In a later paper [32], Nair and coworkers noted that the multifractality of combustion noise tends to disappear as the system transitions to instability and suggested that a decline in the Hurst exponent is a good indicator of this phenomenon. Further studies have demonstrated that measures derived from symbolic time series analysis (STSA) [39] and complex networks [31] are also able to capture the onset of instability.

Machine learning tools have been employed to extract information from pressure signals instead of relying on hand-crafted precursors of instability. Hidden Markov models constructed from the output of STSA

[21] or directly from pressure measurements [30]

have been used to classify the state of combustors. Hachijo et al

[15] 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 [25]

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


show 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


uses the detrended fluctuation analysis (DFA) spectrum of the pressure signal as input to a random forest and finds that this approach compares favorably to precursors from the literature. Recent works have investigated machine learning methods for the design and operation of cryogenic rocket engines

[8, 48, 49, 50].

In this paper, we investigate the question of whether the features of combustion noise are sufficient to reliably predict the occurrence of instabilities in a cryogenic rocket thrust chamber. Certain combustion noise characteristics are believed to be widely unaffected by the geometrical boundary conditions of an engine and thus enable applicability to a wide range of combustion devices with little or no adaption.

Our main contributions are the following:

  • development of a data-driven method to construct early warning signals for thermoacoustic instabilities using nonlinear combustion noise features and support vector machines

  • quantitative evaluation of the constructed early warning signals on data from a representative cryogenic rocket thrust chamber

  • comparison with state-of-the-art early warning indicators

The remainder of the paper is structured as follows: Section II describes the research thrust chamber "D" (BKD) and the experimental data. The basics of recurrence quantification analysis are outlined in section III. Section IV discusses support vector machines. Section V and section VI present the data-driven method for construction of early warning signals and the test case respectively. The results, including the comparison with the performance of other early warning indicators, are shown in section VII. Finally, section VIII provides concluding remarks.

Ii Rocket Thrust Chamber "BKD"

This study is performed with experimental data of the research thrust chamber "D" (BKD), which is operated at the European Research and Technology test bench P8 at the DLR Institute of Space Propulsion in Lampoldshausen. BKD has been designed for heat transfer research with regenerative cooling of liquid hydrogen (LH) at representative conditions of European industrial LOX/H engines, such as Vulcain or Vinci [42, 18]. However, the thrust chamber showed self-excited high-frequency combustion instabilities during the first tests for LOX/H combustion. For that reason, BKD has become a valuable platform for the analysis of chamber acoustic phenomena and the underlying coupling mechanism due to the representative geometry and conditions [13, 10, 11, 12, 18, 1, 2]. In later tests, BKD was also used with a similar setup for the investigation of LOX/LNG combustion, which also showed high-frequency combustion instabilities [27]. However, the current study only focuses on the prediction of thermoacoustic instabilities for the propellant combination LOX/H in BKD.

Figure 1: Experimental thrust chamber BKD [10].

The BKD setup is shown in Fig. 1. It consists of three main components: a multi-element injector head, a cylindrical combustion chamber, and a convergent-divergent nozzle. The injector head has 42 shear coaxial injection elements. The cylindrical combustion chamber has an inner diameter of 80 mm and is 200 mm long. For the instability investigations, the combustion chamber is water-cooled in order to guarantee a sufficient safety margin with respect to increasing thermal loads. The nozzle has a throat diameter of 50 mm, which leads to representative chamber characteristics as a contraction ratio of and a characteristic chamber length of  m.

For the investigation of the stability behavior of this particular device, a large number of different operating conditions have been tested in several test campaigns. An operating condition or load point (LP) is thereby primarily defined by the mean chamber pressure , the propellant mixture ratio (ROF = ) and the propellant injection temperatures ( and ). The chamber pressure has been varied between 50 and 80 bar. For that reason, all tested chamber pressures are close to or above the critical pressure of oxygen, at which LOX injection and combustion is transcritical. ROF between 2 and 7 has been achieved, while most load points are between ROF 4 and 6. During the operation of BKD at the test bench P8, the LOX injection temperature is typically around 110 K. The hydrogen injection temperature depends on the used H storage system of the test bench [14] and is usually around 100 K. If a cryogenic storage tank is used, a lower injection temperature of 45 to 50 K can be achieved. For LOX/H it has been reported that the hydrogen injection temperature has an impact on combustion stability [47, 36, 52]. For that reason, hydrogen temperature ramping tests have been performed in the past [11]. In this case is varied between 45 and 135 K in several ramps.

For the load point of 80 bar ROF 6 the thrust chamber achieves a theoretical thermal power of almost 100 MW and a thrust of about 24 kN. These numbers place BKD at the lower end of small upper stage engines.

The representative conditions in the thrust chamber with temperatures of up to 3600 K and high pressures limit the diagnostics for the instability investigations. A special measurement ring is placed between the injector head and the cylindrical combustion chamber segment, as indicated in Fig. 1. At this location, the temperatures are still moderate due to the injection of cryogenic propellants, which allows the mounted sensors to survive the harsh conditions for several test runs. This measurement ring is extensively equipped with sensors, such as thermocouples or pressure sensors to measure the mean chamber pressure. All sensors are mounted in a common measurement plane, which is positioned 5.5 mm downstream of the injection plane. The main diagnostics for the stability analysis are water-cooled high-frequency piezoelectric pressure sensors. Eight Kistler type 6043A are flush-mounted in the ring with an even circumferential distribution, in order to measure the chamber pressure oscillations . The high-frequency pressure sensors have a measurement range set to  bar and a sampling rate of 100 kHz. An anti-aliasing filter was set at 30 kHz.

Two different types of self-excited high-frequency combustion instabilities have been detected, both of which are driven by injection coupling [13, 10, 18, 1]. A typical BKD test sequence, is shown in Fig. 2. A pressure oscillation spectrogram from inside the combustion chamber is also presented.

Figure 2: BKD test sequence and spectrogram showing the type 1 instability [10].

Stable and unstable operating conditions can be identified in the spectrogram. The strongest self-excited high-frequency combustion instabilities of the first tangential (1T) resonance mode at about 10 kHz were found for the 80 bar, ROF 6 load point between 16 and 25 s with a hydrogen injection temperature of  K. The maximum peak-to-peak amplitudes during unstable combustion reached up to 16-35 bar but showed a large temporal variation. In addition to the pressure oscillation data, Gröning et al. [10, 12] used fibre-optical probes to record fluctuating OH* radiation intensity of individual flames. These optical measurements are not included in the current study but helped to identify the instability driving mechanism in BKD. The OH* fluctuations showed dominant frequencies corresponding to LOX post acoustic resonance frequencies, independent of chamber acoustics [10, 13]. For that reason, Gröning et al. described the coupling mechanism as injection-driven [10]. The flame dynamics are modulated by the LOX post acoustics, and combustion instabilities emerge when the frequency of a chamber mode matches one of the longitudinal modes of the injectors. For the load points of 80 bar ROF 6 this frequency interaction appears for the chamber 1T mode with the second longitudinal eigenmodes of the LOX posts [10]. This mechanism was later confirmed with high-speed imaging of the flame dynamics in the thrust chamber using a 2D optical access window [2]. Recent investigations indicated that an additional hydrodynamic effect in the LOX injectors may also play a role in the excitation of the LOX post eigenmodes and the amplification of the combustion instability [2]. As was shown, the coupling mechanism of this type of instability seems well understood. Furthermore, the stability behavior shows reproducible characteristics. If the load point of 80 bar ROF 6 was approached in the test sequence, the chamber consistently showed an excitation of the 1T mode for the identical thrust chamber configuration [10, 13].

Beside the aforementioned LOX-coupled type instability, a second type with very high amplitudes (up to > 75 % of ) appears under rare circumstances [18, 1]. This type of instability appeared for different operating conditions and injector lengths. Most of the load points, which showed this type of instability, have a cold hydrogen injection temperature of about 45 K. Fig. 3 shows a test sequence and -spectrogram of a run with cold hydrogen injection, which showed this type of instability.

Figure 3: BKD test sequence and -spectrogram showing the type 2 instability [1]

As can be observed, the instability shows higher amplitudes compared to the type 1 instability, sometimes even exceeding the measurement range of the acoustic pressure sensors, and appears for different operating conditions. Furthermore, even for stationary operating conditions the high-amplitude pressure oscillations can appear and disappear spontaneously several times, as can, for example, be observed for the load point of 70 bar ROF 6 during the time window from 26-30 s. Other effects, which can be observed in Fig. 3 are an increase of the hydrogen injection temperature during the instability and an apparent broadening of the acoustic spectrum due to switching of instability type and therefore frequency. Recent investigations of the acoustic field dynamics and the coupling mechanism indicated a three-way coupling with the chamber and both injection systems [1]. At first, the instability seems to be initialized similarly to the type 1 instability, by a coupling of the LOX-posts with the chamber acoustics. However, for unknown reasons, the amplitudes grow significantly higher than that of the type 1 instability. The resulting strong transverse oscillations enhance the mixing of the injection propellants, which leads to short flames [1]. Similar observations have been reported in other experiments [17, 16] and simulations [18, 3, 16, 40, 45]. Due to the shortened flames, the 1T mode frequency increases and allows coupling with the H injectors [1].

The fact that this type of oscillation appears spontaneously during stationary operating conditions indicates that the system is being triggered by a nonlinear mechanism. Previous studies were not able to identify a triggering process although in gas turbines, triggering has been seen to be caused by background noise [53], due to the non-normal behaviour of thermoacoustic systems [23]. For that reason, while being able to identify operating conditions with an increased risk for this type of instability, it is currently impossible to predict the onset of this instability. These characteristics make the prediction and analysis of the type 2 instability more complicated than the type 1 instability.

Iii Recurrence Quantification Analysis

In this section, we review the basics of recurrence quantification analysis (RQA) [28, 51]. RQA was developed to study and compare recurrence plots which are used to visualize the recurrence behavior of the phase space trajectory of dynamical systems. By Takens’ delay embedding theorem [44], we can reconstruct the dynamics of the rocket thrust chamber in an appropriate phase space from a single state variable such as the acoustic pressure [22]. The univariate pressure time series data are converted into a set of delayed vectors:


where and denote an appropriate time delay and embedding dimension respectively. The time delay

can be estimated using the autocorrelation function or average mutual information

[24]. The embedding dimension can be obtained using the false nearest neighbor method or Cao’s method [5]. After the reconstruction of a suitable phase space, the recurrence matrix is computed by the following formula:


where is the Heaviside step function, is a predefined threshold, and is the Euclidean distance between the state points and . Whenever a state point recurs in the predefined threshold, is set to 1, otherwise is set to 0. A recurrence plot is the two-dimensional arrangement of recurrence points, i.e. a visualisation of the recurrence matrix. Different patterns characterize different dynamics of the signal. Several statistical measures are used to quantify the structures present in a recurrence plot [28, 51]. The recurrence rate (RR) measures the density of recurrence points:


where is the number of state vectors in the reconstructed phase space. The determinism (DET) measures the percentage of recurrence points which form diagonal lines of minimal length :



is the probability distribution of diagonal lines having length

. The measure is called determinism because it is related to the predictability of the studied dynamical system. The laminarity (LAM) similarly quantifies the amount of recurrence points which form vertical lines:



is the probability distribution of vertical lines having length

, which have at least a length of . The probability that a diagonal line has exactly length can be estimated from the probability distribution :


Using , one can further calculate the associated Shannon entropy (ENTR):


which reflects the complexity of the deterministic structure in the dynamical system. An additional measure is the ratio of the determinism to the recurrence rate, i.e. RATIO=DET/RR.

Iv Support Vector Machines

In this section, we briefly recall the basics of support vector machines (SVMs) [46, 4]

. SVMs are machine learning models with associated (supervised) learning algorithms, which are used for classification and regression tasks. In binary classification, one is given a training dataset of

points of the form


with -dimensional patterns and associated class labels , and tries to estimate a function such that will correctly classify new examples , i.e. for points which were generated from the same underlying probability distribution as the training data. SVMs assume that the decision function

is based on the class of hyperplanes


where and . The associated function is of the form


The goal is to find the hyperplane with the maximal margin of separation between the two classes. In general, the larger the margin, the lower the generalization error of the classifier. If the training data is linearly separable, this leads to the following optimization problem:


for . The above is an optimization problem with a convex quadratic objective and only linear constraints. It can be solved using quadratic programming techniques. Recent methods for finding the SVM classifier also include sub-gradient descent and coordinate descent algorithms. The optimal hyperplane is completely determined by those that lie nearest to it. These data points are called support vectors. It follows that is given by a linear combination of the support vectors,


and the decision function can be written as


To extend SVMs to cases in which the data are not linearly separable, one introduces slack variables to allow certain constraints to be violated. The new optimization problem can be expressed as


and , for all . Nonlinear classifiers can be constructed by applying the kernel trick below. The resulting algorithm is formally similar, except that every dot product is replaced by a nonlinear kernel function


where is a nonlinear map from into a feature space , . Although the algorithm fits a hyperplane in the transformed feature space, the decision function may be nonlinear in the original input space,


The hyperparameters consist of the soft margin constant,

, and parameters on which the kernel function may depend, e.g. the width of a Gaussian kernel.

V Early Warning Signal Construction

In this section, we present the proposed construction of a suitable early warning signal. Clearly, a suitable warning signal could be constructed using a combustion noise feature that increases or decreases monotonically before an instability and a universal threshold, which separates stable data points depending on whether they belong to the onset of instability or not. The comparison of the current value with the threshold could be used as an early warning signal. In general, this will not be possible and the classification using a simple threshold value for a certain feature will not correspond to a division into data points that belong to the onset of instability or not. Naively, one could think that the magnitude of the amplitude spectrum of the pressure signal at the dominant frequency constitutes a suitable measure by setting an appropriate threshold. A closer look shows that this measure is not well suited, because of the dependency on operating conditions, propellant combinations, and combustion chamber geometries [24]. Furthermore, the root mean square (RMS) of the pressure signal does not increase in a monotonic way as the system approaches a combustion instability [24]. Thus, for such features, it is not possible to define universal thresholds to detect the onset of thermoacoustic combustion instability.

To overcome the shortcomings of conventional measures, different measures from nonlinear time series analysis have been introduced. Certain nonlinear combustion noise features are more independent of operating conditions and the exact combustion chamber geometries. RQA measures quantify recurrence behaviors, which are different for stable combustion and combustion instability [22]

. The pressure signals in the unstable (limit-cycle) regime possess a deterministic periodic nature, while the stable regime is distinguished by a noisy or chaotic nature. Furthermore, the transient phase is characterized by a strong change in these measures. Because of this, it makes sense not only to look at the current values, which describe the present recurrence behavior, but also to investigate their variation. Trend estimation techniques can be used to determine whether time series data exhibit an increasing or decreasing trend. The linear fit for different sample windows is one of the simplest methods. It is expected that the consideration of trend behavior should increase predictability. A combination of several combustion noise features reduces the influence of outliers and usually makes forecasts more stable. To find the optimal combination and decision criterion respectively, one can use data-driven machine learning algorithms.

Thus, we apply the following procedure:

  • calculate RQA measures and estimate their trends using suitable sample windows

  • train an SVM to learn the associated magnitudes and trends in the data that belong to different regimes

  • use the class prediction of the trained SVM as an online early warning signal

Compared to the use of only one combustion noise feature, the proposed method is expected to generate fewer false alarms and has the potential to work with other combustion chamber geometries due to the use of more universal combustion noise features. The central challenge is that the machine learning model may overfit during training by memorizing properties of the training data that do not work well on unseen data [4]. The capability to perform well on unseen input data is called generalization and can be estimated by using only a fraction of the available data for training. The remaining portion is used to select the optimal hyperparameter combination.

Vi Test Case

For the evaluation, we use a data set of 10 typical BKD tests. Table 1 summarizes the main characteristics of each run.

Run Propellant state111GH (gaseous H2) and LH (liquid H2) refer to different fuel interfaces of the test bench. At the LH interface the hydrogen is stored as a cryogenic fluid, while at the GH interface it is stored in a high-pressure tank under ambient conditions and a heat exchanger is to used to further cool it down. For both cases, hydrogen is injected as a supercritical fluid, but the injection temperatures vary. The LH interface leads to an injection temperature around 45 K. The GH interface results in an injection temperature around 100 K. Post length222Most experiments used a LOX-post with a length of 68 mm. Additional tests with modified lengths, e.g. 64 mm, were also carried out. Instability type
1 GH 68 1
2 GH 68 1
3 GH 68 1
4 GH 68 1
5 GH 68 1
6 LH 68 1,2
7 GH 64 2,1
8 GH 68 2,1,1
9 LH 68 1,2
10 LH 68 1,2
Table 1: Overview of used BKD test runs. Fig. 11 and Fig. 12 in the appendix show the growth of the normalized peak-to-peak amplitudes of the pressure oscillations and the associated test sequences.

In total there are 16 instabilities present in the set of time series. The threshold of a data point to be classified as type 1 instability is set to a peak-to-peak amplitude of 6.25 % with respect to the mean chamber pressure. This threshold marks the transition to limit cycle oscillations in BKD. A threshold of 5 % would be violated by single spikes, which are not part of a lasting instability. In terms of type 2 instabilities, a threshold of 20.0 % is chosen due to their large amplitudes. Fig. 11 and Fig. 12 in the appendix show the growth of the normalized peak-to-peak amplitudes and the associated test sequences. The training and validation data set is given by the BKD runs . The test data set is given by the remaining BKD runs .

First, for each time series, i.e. pressure signal of a BKD run, we calculate the points in time when an instability begins and when an instability ends. We define the start of an unstable combustion process by the fact that the peak-to-peak amplitude increases above 6.25 % of the mean chamber pressure. An instability ends when the value drops below 6.25 % of the mean chamber pressure again and stays there for at least 500 ms. The second condition causes a fast sequence of stable and unstable phases to be classified as unstable. This is desired because we intend to predict the first occurrence of strong pressure fluctuations.

Second, the following RQA measures are calculated using a sliding window approach: RR, DET, LAM, ENTR, and RATIO. The calculation of the quantities for time uses the values of the dynamic pressure signal from the interval . For each measure, we also estimate the time-dependent slope of the trend using a linear fit and windows of the form . Further details of the calculation are described in the appendix. In total, we obtain ten derived combustion noise features.

To characterize the transient regime, the data points which belong to the stable phase are divided into two sets. The first set contains all stable data points whose points in time are at least 200 ms away from the occurrence of instability, while the second set includes all stable data points which are not in the first set. An SVM is used to solve this binary classification problem. We use a random search to find the best hyperparameters. To compare different hyperparameter combinations, we use cross-validation to assess the performance of the SVM measured by the F-score (harmonic mean of precision and recall). Data from one run are used to evaluate the prediction accuracy, while data from the other runs are used to train the SVM for a given hyperparameter combination. This procedure is repeated, such that data from every run are used for evaluation once. Finally, the performance is averaged. The hyperparameters which belong to the best average performance are used to train the final SVM using the complete training and validation data set. Finally, the prediction of this SVM is investigated using the test data set.

Vii Results

This section presents the results of the presented data-driven early warning signal method for the BKD test case. The performance of the classifier reflects the capability to successfully predict if the system will develop a thermoacoustic instability within the next 200 ms using dynamic pressure sensor data. In other words, we quantify if the model can detect a divergence from the stable regime. Online prediction with such a time window provides a realistic lead-time for taking appropriate control actions.

For a classification problem with imbalanced classes, i.e. where the number of samples in the classes differs greatly as in our case, it is crucial to use suitable performance measures. The true-positive rate (TPR) is given by the number of true positives, i.e. the number of data points correctly labeled as belonging to the positive class (transient regime in our case), divided by the total number of data points that belong to the positive class. The false-positive rate (FPR) is given by the number of false positives, i.e. the number of data points incorrectly labeled as belonging to the positive class, divided by the total number of data points that belong to the negative class (far from instability in our case). A high FPR is equivalent to a large number of false alarms, i.e. a warning signal wrongly predicting a thermoacoustic instability. Thus, for suitable early warning signals, it is extremely important to exhibit a small FPR. Table 2 shows the TPRs associated with the trained SVM for different FPRs. For an FPR of 0.5 % the TPR is 49 %. The TPR increases to 61 % and 80 % for an FPR of 1 % and 2 % respectively. A receiver operating characteristic (ROC) curve, can be used to illustrate the diagnostic ability of a binary classifier as its discrimination threshold is varied. It is created by plotting the TPR against the FPR at various threshold settings. Fig. 4 displays the ROC curve for the SVM. One can see that the TPR remains limited to about 80 % even for large FPRs.

Figure 4: Receiver operating characteristic (ROC) curve of the test data set for the trained SVM classifier. Points above the diagonal represent good classification results (better than random).

By looking at the decision function of the SVM for all instabilities in the test set, the responsible reasons become apparent. The decision function and the threshold associated with an FPR of 1 % is shown in Fig. 5 and Fig. 6.

TPR(FPR=0.5 %) TPR(FPR=1 %) TPR(FPR=2 %)
SVM 49 % 61 % 80 %
Table 2: True-positive rate (TPR) of the SVM classifier for different false-positive rates (FPRs). The FPRs are realized by varying the decision threshold.
Figure 5: Predicted decision function values and the threshold corresponding to an FPR of 1 % near the onset of combustion instability for the first part of the test data set. For evaluation of the early warning signal, the normalized peak-to-peak amplitude of the pressure oscillations is also shown. The start of combustion instability is defined by the condition that the normalized peak-to-peak amplitude increases above 6.25 % (type 1) and 20.0 % (type 2) respectively. An instability ends when the value drops below 6.25 % again and stays there for at least 500 ms. The area of instability is marked red.
Figure 6: Predicted decision function values and the threshold corresponding to an FPR of 1 % near the onset of combustion instability for the second part of the test data set. For evaluation of the early warning signal, the normalized peak-to-peak amplitude of the pressure oscillations is also shown. The start of combustion instability is defined by the condition that the normalized peak-to-peak amplitude increases above 6.25 % (type 1) and 20.0 % (type 2) respectively. An instability ends when the value drops below 6.25 % again and stays there for at least 500 ms. The area of instability is marked red.

The instability in run 3 is predicted around 200 ms before occurrence. The type 2 instability in run 7 is predicted within 200 ms, but the value of the decision function exceeds the threshold 400 ms before the instability begins. Thus, there are false positives. The transition to the type 1 instability is detected as intended. Nevertheless, the performance for this run is very convincing, because run 7 was carried out with a modified injector geometry and there are no time series for this setup in the training data. The fact that the prediction still performs so well indicates that the early warning signal could also work with modified injection systems. The performance in run 8 is not convincing. Only the second instability, which is of type 1, is satisfactorily predicted. At the first instability (type 2) the warning signal flickers but goes out before the instability. The third instability of type 1 is detected too late. The performance at the last run in the test set is again very convincing and the early warning signal works as intended. This is astonishing because the time series belongs to an experiment with a reduced injection temperature of the hydrogen.

All in all, in 6 of 8 instabilities (type 1 and type 2) the constructed early warning signal works in a satisfying way. This leads to a TPR of about 60 %. The TPR can be improved by reducing the threshold, but this also increases the number of false alarms. To predict the instabilities that are present at run 8 in time, the threshold must be significantly reduced. This leads to FPRs that are no longer acceptable.

Figure 7: Permutation feature importance plot. Feature importances are measured through the decrease of the F-score on the test data set.

To estimate the relative importance of input features, one can use the permutation feature importance method. Permutation feature importance measures the decrease in the prediction accuracy of the model after permuting the feature. A feature is more important if shuffling its values decreases the performance to a greater extent. Fig. 7 shows the feature importances measured through the features’ influence on the F-score. The most important input features are given by the recurrence rate RR and the Shannon entropy ENTR.

For comparison, we also calculate the TPRs and FPRs, which result from the use of a single measure. By varying the threshold, the TPR can be determined for different FPRs. The results are shown in Table 3. For an FPR of 1 % the TPRs are less than or equal to 40 %. Our proposed method leads to a TPR of 61 %.

Measure TPR(FPR=0.5 %) TPR(FPR=1 %) TPR(FPR=2 %)
RR 6 % 32 % 59 %
DET 14 % 35 % 46 %
LAM 2 % 40 % 71 %
ENTR 32 % 36 % 38 %
RATIO 2 % 33 % 62 %
Table 3: True-positive rates (TPRs) associated with single RQA measures for different false-positive rates (FPRs). The FPRs are realized by varying the decision thresholds.

Another measure that can be calculated from the combustion noise is the Hurst exponent . There is a loss of multifractality in combustion noise as combustors progress towards combustion instability, which is reflected in a decline of [32, 22]. We compute by calculating the slope of the detrended fluctuations for logarithmically spaced window sizes ranging from 50 000 to 70 000 data points (corresponding to 5 to 7 cycles of the 10 kHz instability). In the paper of Nair et al [32], it is suggested to use 2 to 4 cycles, but the slope is changing a lot in this range. Thus, we increase the number of cycles to get the asymptotic slope of the fluctuations. Although drops to a very low value when the instability is present, the performance as an early warning signal is not convincing. Table 4 shows the prediction accuracy using an optimal threshold for the test data set. When using 2 to 4 cycles for the computation of , the prediction accuracy is even lower. The main reason for this is that fluctuates strongly and therefore generates many false alarms. As an example, Fig. 8 displays the time development of for run 7.

TPR(FPR=0.5 %) TPR(FPR=1 %) TPR(FPR=2 %)
11 % 24 % 27 %
Table 4: True-positive rate (TPRs) associated with the Hurst exponent for different false-positive rates (FPRs). The FPRs are realized by varying the decision thresholds.
Figure 8: Hurst exponent for run 7.

Viii Conclusion and Outlook

In this paper, we presented the data-driven construction of a warning signal for the early detection of thermoacoustic instabilities in rocket thrust chambers. Using time-series data from experimental tests performed with the cryogenic rocket thrust chamber BKD, we calculated time-dependent RQA measures and their variation. Then we applied machine learning to obtain the characteristic properties of the transition phase. After the training and tuning phase, we evaluated the method in a quantitative way and compared it with other early warning indicators.

The constructed early warning signal achieved the best performance and was able to predict 6 of 8 thermoacoustic instabilities (type 1 and type 2) 200 ms in advance with a low probability of false alarms. Furthermore, the results indicate that there is a good transferability to other injection systems. There might be a good transferability to other combustion chamber geometries. With this approach, high-frequency combustion instabilities can be successfully predicted under representative conditions. This is an important step towards the active control of combustion instabilities in the future.

There are essential questions that should be investigated in the future. Besides the systematic examination of the generalization ability, the influence of different noise levels should be examined. Since many research thrust chambers like the BKD are equipped with diverse sensors, methods should also be investigated, which fuse data from multiple sensors.

Information contained in the pressure signal is lost in the calculation of the RQA measures. Therefore machine learning methods should be investigated which derive suitable features directly from the data, e.g. deep neural networks.

Machine learning techniques can perform poorly when given inputs dissimilar to those seen during training, so-called out-of-distribution data. This can be problematic for safety-critical systems like rocket engines, particularly since our algorithms are trained on data sets whose size is necessarily limited by the resource-intensive nature of the experiments. We are therefore exploring Bayesian neural networks as a tool to quantify predictive uncertainties and make our predictions robust to overconfident extrapolations.

The authors would like to thank the crew of the P8 test bench and also Stefan Gröning for preparation and performing the test runs. Furthermore, it is a pleasure to thank Kai Dresia and Christoph Räth for many useful discussions.

Appendix A Calculation Details

The basis of all calculations is the signal of a single dynamic pressure sensor. A high-pass filter is used to remove non-acoustic contributions. For reconstructing the phase space that is traced by the pressure oscillation, we estimate the optimum time delay and embedding dimension using the autocorrelation function [35] and Cao’s method [5] respectively. The optimum time delay corresponds to the first zero crossing of the autocorrelation function and is given by 0.02 ms. We chose a constant embedding dimension of 15 for all time series. For details related to the method of Cao to identify the minimum embedding dimension, the reader is referred to the paper of Cao [5]. The time-dependent recurrence plots use the signal values from the last 200 ms in a sliding window manner. The threshold is set to 3 distance units after rescaling the input. Fig. 9 and Fig. 10 show two exemplary recurrence plots. The pattern that is visible in Fig. 9 is typical for a stable combustion process. The pattern that is shown in Fig. 10 characterizes a transition to instability.

Figure 9: Recurrence plot for the typical dynamics of stable operation.
Figure 10: Recurrence plot for the typical dynamics of a transition to instability.

By using recurrence plots analogous to the shown ones, RQA measures are calculated using the pyunicorn Python package [7]

. For trend estimation, we employ values of the last 100 ms in a sliding window manner and a simple linear fit. Then, all input features are standardized by removing the mean and scaling to unit variance. We use the standard scikit-learn

[37] SVM implementation for the binary classification and adjust the class weights inversely proportional to class frequencies in the input data. The optimal hyperparameters turn out to be a soft margin constant

and a radial basis function kernel coefficient


Appendix B BKD Runs

Fig. 11 and Fig. 12 show the normalized peak-to-peak amplitudes and the associated test sequences for all BKD runs used in the current study.

Figure 11: Normalized peak-to-peak amplitudes of the pressure oscillations and the associated test sequences for the BKD runs used in the current study. The area of instability is marked red.
Figure 12: Normalized peak-to-peak amplitudes of the pressure oscillations and the associated test sequences for the BKD runs used in the current study. The area of instability is marked red.


  • [1] W. Armbruster, J. S. Hardi, D. Suslov, and M. Oschwald (2018) Experimental investigation of self-excited combustion instabilities with injection coupling in a cryogenic rocket combustor. Acta Astronautica 151, pp. 655–667. External Links: Document Cited by: Figure 3, §II, §II, §II, §II.
  • [2] W. Armbruster, J. S. Hardi, D. Suslov, and M. Oschwald (2019) Injector-driven flame dynamics in a high-pressure multi-element oxygen-hydrogen rocket thrust chamber. Journal of Propulsion and Power 35, pp. 632–644. External Links: Document Cited by: §II, §II.
  • [3] S. K. Beinke (2017) Anaylses of flame response to acoustic forcing in a rocket combustor. Ph.D. Thesis, School of Mechanical Engineering, The University of Adelaide, Australia. External Links: Document Cited by: §II.
  • [4] C. M. Bishop (2006) Pattern recognition and machine learning. Springer. Cited by: §IV, §V.
  • [5] L. Cao (1997) Practical method for determining the minimum embedding dimension of a scalar time series. Physica D: Nonlinear Phenomena. External Links: Document Cited by: Appendix A, §III.
  • [6] S. Colas, S. L. Gonidec, P. Saunois, M. Ganet, A. Remy, and V. Leboeuf (2019) A point of view about the control of a reusable engine cluster. In Proceedings of the 8th European Conference for Aeronautics and Space Sciences, Madrid, Spain. External Links: Document Cited by: §I.
  • [7] J. F. Donges, J. Heitzig, B. Beronov, M. Wiedermann, J. Runge, Q. Y. Feng, L. Tupikina, V. Stolbova, R. V. Donner, N. Marwan, H. A. Dijkstra, and J. Kurths (2015) Unified functional network and nonlinear time series analysis for complex systems science: The pyunicorn package. Chaos. External Links: Document Cited by: Appendix A.
  • [8] K. Dresia, G. Waxenegger-Wilfing, J. Riccius, J. Deeken, and M. Oschwald (2019) Numerically Efficient Fatigue Life Prediction of Rocket Combustion Chambers using Artificial Neural Networks. In Proceedings of the 8th European Conference for Aeronautics and Space Sciences, Madrid, Spain. External Links: Document Cited by: §I.
  • [9] H. Gotoda, H. Nikimoto, T. Miyano, and S. Tachibana (2011) Dynamic properties of combustion instability in a lean premixed gas-turbine combustor. Chaos. External Links: Document, ISSN 10541500 Cited by: §I.
  • [10] S. Gröning, J. S. Hardi, D. Suslov, and M. Oschwald (2016) Injector-driven combustion instabilities in a Hydrogen/Oxygen rocket combustor. Journal of Propulsion and Power 32 (3), pp. 560–573. External Links: Document Cited by: Figure 1, Figure 2, §II, §II, §II.
  • [11] S. Gröning, J. S. Hardi, D. Suslov, and M. Oschwald (2017) Influence of hydrogen temperature on the stability of a rocket engine combustor operated with hydrogen and oxygen. CEAS Space Journal 9 (1), pp. 59–76. External Links: Document Cited by: §II, §II.
  • [12] S. Gröning, J. S. Hardi, D. Suslov, and M. Oschwald (2019) Measuring the phase between fluctuating pressure and flame radiation intensity in a cylindrical combustion chamber. Progress in Propulsion Physics 11, pp. 425–446. External Links: Document Cited by: §II, §II.
  • [13] S. Gröning (2017) Untersuchung selbsterregter verbrennungsinstabilitäten in einer raketenbrennkammer. Ph.D. Thesis, RWTH Aachen. Cited by: §II, §II, §II.
  • [14] A. Haberzettl, D. Suslov, G. Brümmer, A. Frank, and M. Oschwald (2016) European technology test facility for high pressure combustion P8. 20 years reserach and technology development. In Deutscher Luft- und Raumfahrtkongress (DLRK) 2016, Braunschweig, Germany. Cited by: §II.
  • [15] T. Hachijo, S. Masuda, T. Kurosaka, and H. Gotoda (2019) Early detection of thermoacoustic combustion oscillations using a methodology combining statistical complexity and machine learning. Chaos: An Interdisciplinary Journal of Nonlinear Science 29 (10), pp. 103123. External Links: Document, Link, https://doi.org/10.1063/1.5120815 Cited by: §I.
  • [16] L. Hakim, A. Ruiz, T. Schmitt, M. Boileau, G. Staffelbach, S. Ducruix, B. Cuenot, and S. Candel (2015) Large eddy simulations of multiple transcritical coaxial flames submitted to a high-frequency transverse acoustic modulation. Proceedings of the Combustion Institute 35 (2), pp. 1461–1468. External Links: Document Cited by: §II.
  • [17] J. S. Hardi, S. K. Beinke, M. Oschwald, and B. B. Dally (2014) Coupling of cryogenic oxygen"=hydrogen flames to longitudinal and transverse acoustic instabilities. Journal of Propulsion and Power 30 (4), pp. 991–1004. External Links: Document Cited by: §II.
  • [18] J. S. Hardi, T. Traudt, C. Bombardieri, M. Börner, S. K. Beinke, W. Armbruster, P. N. Blanco, F. Tonti, D. Suslov, B. Dally, and M. Oschwald (2018) Combustion dynamics in cryogenic rocket engines: research programme at dlr lampoldshausen. Acta Astronautica 147, pp. 251–258. External Links: Document Cited by: §II, §II, §II, §II.
  • [19] D. Harrje and F. Reardon (Eds.) (1972) Liquid propellant rocket combustion instability. NASA SP-194. Cited by: §I, §I.
  • [20] S. D. Heister, W. E. Anderson, T. L. Pourpoint, and R. J. Cassady (2019) Rocket propulsion. 1 edition, Cambridge University Press. Cited by: §I.
  • [21] D. K. Jha, N. Virani, and A. Ray (2018) Markov modeling of time series via spectral analysis for detection of combustion instabilities. In Handbook of Dynamic Data Driven Applications Systems, pp. 123–138. Cited by: §I.
  • [22] M. P. Juniper and R. I. Sujith (2018) Sensitivity and nonlinearity of thermoacoustic oscillations. Annual Review of Fluid Mechanics. External Links: Document Cited by: §I, §III, §V, §VII.
  • [23] M. P. Juniper (2011) Triggering in the horizontal rijke tube: non-normality, transient growth and bypass transition. Journal of Fluid Mechanics. External Links: Document Cited by: §II.
  • [24] P. Kasthuri, I. Pavithran, S. A. Pawar, R. I. Sujith, R. Gejji, and W. Anderson (2019) Dynamical systems approach to study thermoacoustic transitions in a liquid rocket combustor. Chaos. External Links: Document Cited by: §III, §V.
  • [25] T. Kobayashi, S. Murayama, T. Hachijo, and H. Gotoda (2019-06) Early detection of thermoacoustic combustion instability using a methodology combining complex networks and machine learning. Phys. Rev. Applied 11, pp. 064034. External Links: Document, Link Cited by: §I.
  • [26] T. Lieuwen (2005) Online Combustor Stability Margin Assessment Using Dynamic Pressure Data. Journal of Engineering for Gas Turbines and Power. External Links: Document, ISSN 07424795 Cited by: §I.
  • [27] J. Martin, W. Armbruster, J. Hardi, D. Suslov, and M. Oschwald (2020) Experimental investigation of self-excited combustion instabilities in a lox/lng rocket combustor. In AIAA Scitech 2020 Forum, AIAA Paper 2020-0423, Orlando, Florida. External Links: Document Cited by: §II.
  • [28] N. Marwan (2003) Encounters with neighbours: current developments of concepts based on recurrence plots and their applications. Norbert Marwan. Cited by: §III.
  • [29] M. McCartney, T. Indlekofer, and W. Polifke (2020) Online prediction of combustion instabilities using machine learning. Journal of Engineering for Gas Turbines and Power. Cited by: §I.
  • [30] S. Mondal, N. F. Ghalyan, A. Ray, and A. Mukhopadhyay (2018) Early detection of thermoacoustic instabilities using hidden markov models. Combustion Science and Technology. Cited by: §I.
  • [31] M. Murugesan and R. Sujith (2016) Detecting the onset of an impending thermoacoustic instability using complex networks. Journal of Propulsion and Power 32 (3), pp. 707–712. Cited by: §I.
  • [32] V. Nair and R. I. Sujith (2014-05) Multifractality in combustion noise: predicting an impending combustion instability. Journal of Fluid Mechanics 747, pp. 635–655. External Links: Document, ISSN 0022-1120, Link Cited by: §I, §VII.
  • [33] V. Nair, G. Thampi, S. Karuppusamy, S. Gopalan, and R. I. Sujith (2013) Loss of chaos in combustion noise as a precursor of impending combustion instability. International Journal of Spray and Combustion Dynamics 5 (4), pp. 273–290. External Links: Document, Link, https://doi.org/10.1260/1756-8277.5.4.273 Cited by: §I.
  • [34] V. Nair, G. Thampi, and R. Sujith (2014) Intermittency route to thermoacoustic instability in turbulent combustors. Journal of Fluid Mechanics 756, pp. 470. Cited by: §I.
  • [35] A. H. Nayfeh and B. Balachandran (2008) Applied nonlinear dynamics: analytical, computational, and experimental methods. John Wiley & Sons. Cited by: Appendix A.
  • [36] Y. Nunome, T. Onodera, M. Sasaki, T. Tomita, K. Kobayashi, and Y. Daimon (2011) Combustion instability phenomena observed during cryogenic hydrogen injection temperature ramping tests for single coaxial injector elements. In 47th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit 2011, AIAA Paper 2011-6027, San Diego, California. External Links: Document Cited by: §II.
  • [37] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay (2011) Scikit-learn: machine learning in Python. Journal of Machine Learning Research 12, pp. 2825–2830. Cited by: Appendix A.
  • [38] J. W. S. Rayleigh (1878) The explanation of certain acoustical phenomena. Nature 18 (455), pp. 319–321. Cited by: §I.
  • [39] S. Sarkar, S. R. Chakravarthy, V. Ramanan, and A. Ray (2016) Dynamic data-driven prediction of instability in a swirl-stabilized combustor. International Journal of Spray and Combustion Dynamics 8 (4), pp. 235–253. Cited by: §I.
  • [40] T. Schmitt, G. Staffelbach, S. Ducruic, S. Gröning, J. Hardi, and M. Oschwald (2017) Large-eddy simulations of a sub-scale liquid rocket combustor: influence of fuel injection temperature on thermo-acoustic stability. In 7th European Conference for Aeronautics and Aerospace Sciences (EUCASS) 2017, Milano, Italy. External Links: Document Cited by: §II.
  • [41] M. Schulze and T. Sattelmayer (2017) Linear stability assessment of a cryogenic rocket engine. International Journal of Spray and Combustion Dynamics 9 (4), pp. 277–298. External Links: Document Cited by: §I.
  • [42] J. Sender, D. I. Suslov, J. Deeken, S. Gröning, and M. Oschwald (2016) "L42" technology demonstrator: operational experience. In Space Propulsion 2016, Rome, Italy. Cited by: §II.
  • [43] U. Sengupta, C. Rasmussen, and M. Juniper (2020) Bayesian machine learning for the prognosis of combustion instabilities from noise. Journal of Engineering for Gas Turbines and Power. Cited by: §I.
  • [44] F. Takens (1980) Detecting strange attractors in turbulence. In Dynamical Systems and Turbulence, Warwick 1980, Warwick, United Kingdom. Cited by: §III.
  • [45] A. Urbano, L. Selle, G. Staffelbach, B. Cuenot, T. Schmitt, S. Ducruix, and S. Candel (2016) Exploration of combustion instability triggering using large eddy simulation of a multiple injector liquid rocket engine. Combustion and Flame 169, pp. 129–140. External Links: Document Cited by: §II.
  • [46] V. N. Vapnik (1999)

    An overview of statistical learning theory

    IEEE Transactions on Neural Networks 10 (5), pp. 988–999. External Links: Document Cited by: §IV.
  • [47] J. P. Wanhainen, H. C. Parish, and E. W. Conrad (1966) Effect of propellant injection velocity on screech in 20000-pound hydrogen-oxygen rocket engine. Technical report NASA TN D-3373, National Aeronautics and Space Administration, National Aeronautics and Space Administration, Cleveland, Ohio (USA). Lewis Research Center, Washington D.C.. Cited by: §II.
  • [48] G. Waxenegger-Wilfing, K. Dresia, J. C. Deeken, and M. Oschwald (2020-04) Heat Transfer Prediction for Methane in Regenerative Cooling Channels with Neural Networks. Journal of Thermophysics and Heat Transfer 34 (2), pp. 347–357. External Links: ISSN 1533-6808, Document Cited by: §I.
  • [49] G. Waxenegger-Wilfing, K. Dresia, J. C. Deeken, and M. Oschwald (2020)

    A Reinforcement Learning Approach for Transient Control of Liquid Rocket Engines

    arXiv:2006.11108 [cs, eess, math, stat]. Cited by: §I.
  • [50] G. Waxenegger-Wilfing, K. Dresia, M. Oschwald, and K. Schilling (2020) Hardware-In-The-Loop Tests of Complex Control Software for Rocket Propulsion Systems. In Proceedings of the 71st International Astronautical Congress, Virtual. Cited by: §I.
  • [51] C. Webber and N. Marwan (2015) Recurrence quantification analysis. Springer. Cited by: §III.
  • [52] V. Yang and W. Anderson (Eds.) (1995) Liquid rocket engine combustion instability. American Institute of Aeronautics and Astronautics, Washington, DC. Cited by: §I, §II.
  • [53] B. T. Zinn and T. C. Lieuwen (2005) Combustion instabilities: basic concepts. Combustion Instabilities in Gas Turbine Engines. Cited by: §II.