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

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.

High-amplitude pressure oscillations known as combustion instability are an issue for all chemical propulsion systems. Combustion instabilities are particularly problematic for rocket thrust chambers because of their high energy release rates. Especially thermoacoustic oscillations are a major hazard, but difficult to predict. An important question is whether features of combustion noise can be used to construct reliable early warning signals for representative rocket thrust chambers. Among other things, instability precursors are needed for active combustion control systems. In this study, we use recurrence quantification analysis and well-known machine learning models to construct an early warning signal that should be able to predict high-amplitude instabilities 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. We describe the test case in detail and also evaluate other measures like the Hurst exponent.

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 1 . 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 a) Electronic mail: guenther.waxenegger@dlr.de b) Also at Institute of Jet Propulsion and Turbomachinery, RWTH Aachen University, 52062 Aachen, Germany 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 2 . Combustion instability refers to highamplitude 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 3 . Those instabilities can be further sub-categorized into high-frequency (screeching) and lowfrequency instabilities 2,4 . 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 5 . 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 6 .
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 2 ) or on the other hand through active control. In the case of active con-trol, 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 7 . The first step in this direction was taken by Lieuwen 8 , 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 9 . Similarly, Gotoda and coworkers have used the Wayland test for nonlinear determinism 10 to assess the stability margin of a combustor. Nair et al 11 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 12 , 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) 13 and complex networks 14 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 handcrafted precursors of instability. Hidden Markov models constructed from the output of STSA 15 or directly from pressure measurements 16 have been used to classify the state of combustors. Hachijo et al 17 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 18 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 19 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 20 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 [21][22][23][24] .
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 2 ) at representative conditions of European industrial LOX/H 2 engines, such as Vulcain or Vinci 25,26 . However, the thrust chamber showed self-excited high-frequency combustion instabilities during the first tests for LOX/H 2 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 [26][27][28][29][30][31][32] . 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 33 . However, the current study only focuses on the prediction of thermoacoustic instabilities for the propellant combination LOX/H 2 in BKD.
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 ε c = 2.56 and a characteristic chamber length of L * = 0.64 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 p cc , the propellant mixture ratio (ROF =ṁ O2 /ṁ H2 ) and the propellant injection temperatures (T O2 and T H2 ). 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 T O2 is typically around 110 K. The hydrogen injection temperature T H2 depends on the used H 2 storage system of the test bench 34 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 2 it has been reported that the hydrogen injection temperature has an impact on combustion stability 4,35,36 . For that reason, hydrogen temperature ramping tests have been performed in the past 29 . In this case T H2 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 p . The highfrequency pressure sensors have a measurement range set to ±30 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 [26][27][28]31 . A typical BKD test sequence, is shown in Fig. 2. A pressure oscillation spectrogram from inside the combustion chamber is also presented.
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 T H2 ≈ 100 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. 28,30 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 27,28 . For that reason, Gröning et al. described the coupling mechanism as injection-driven 28 . 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 p cc 80 bar ROF 6 this frequency interaction appears for the chamber 1T mode with the second longitudinal eigenmodes of the LOX posts 28 . This mechanism was later confirmed with high-speed imaging of the flame dynamics in the thrust chamber using a 2D optical access window 32 . 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 32 . 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 27,28 .
Beside the aforementioned LOX-coupled type instability, a second type with very high amplitudes (up to > 75 % of p cc ) appears under rare circumstances 26,31 . 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 p -spectrogram of a run with cold hydrogen injection, which showed this type of instability.
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 31 . 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 31 . Similar observations have been reported in other experiments 37,38 and simulations 26,[38][39][40][41] . Due to the shortened flames, the 1T mode frequency increases and allows coupling with the H 2 injectors 31 .
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 42 , due to the non-normal behaviour of thermoacoustic systems 43 . 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) 44,45 . 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 46 , 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 p 7 . The univariate pressure time series data are converted into a set of delayed vectors: where τ and d denote an appropriate time delay and embedding dimension respectively. The time delay τ can be esti-mated using the autocorrelation function or average mutual information 47 . The embedding dimension d can be obtained using the false nearest neighbor method or Cao's method 48 . 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 x(i) − x( j) is the Euclidean distance between the state points x(i) and x( j). Whenever a state point recurs in the predefined threshold, R i j is set to 1, otherwise R i j 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 44,45 . The recurrence rate (RR) measures the density of recurrence points: where N = n − (d − 1)τ 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 l min : where P(l) is the probability distribution of diagonal lines having length l. 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: where P(v) is the probability distribution of vertical lines having length v, which have at least a length of v min . The probability p(l) that a diagonal line has exactly length l can be estimated from the probability distribution P(l): Using p(l), 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) 49,50 . 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 n points of the form with d-dimensional patterns x i and associated class labels y i ∈ {±1}, and tries to estimate a function f such that f will correctly classify new examples (x, y), i.e. f (x) = y for points which were generated from the same underlying probability distribution as the training data. SVMs assume that the decision function f is based on the class of hyperplanes where w ∈ R d and b ∈ R. The associated function f 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 i = 1, . . . , n. 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 subgradient descent and coordinate descent algorithms. The optimal hyperplane is completely determined by those x i that lie nearest to it. These data points are called support vectors. It follows that w 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 ξ i to allow certain constraints to be violated. The new optimization problem can be expressed as 14) and ξ i ≤ 0, for all i. 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 R d into a feature space F, Φ : R d → F. 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, C, 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 47 . Furthermore, the root mean square (RMS) of the pressure signal does not increase in a monotonic way as the system approaches a combustion instability 47 . 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 7 . 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  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. 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 50 . 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 I summarizes the main characteristics of each run.
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 [1,2,4,5,6,9]. The test data set is given by the remaining BKD runs [3,7,8,10].
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 t uses the values of the dynamic pressure signal from the interval [t − 200 ms, t]. For each measure, we also estimate the timedependent slope of the trend using a linear fit and windows of the form [t − 100 ms, t]. 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 datadriven 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 II 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.
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.
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 8 FIG. 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. 9 FIG. 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. 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.
To estimate the relative importance of input features, one can use the permutation feature importance method. Permu-FIG. 7. Permutation feature importance plot. Feature importances are measured through the decrease of the F-score on the test data set. tation 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 III. For an FPR of 1 % the TPRs are less than or equal to 40 %. Our proposed method leads to a TPR of 61 %.
Another measure that can be calculated from the combus- tion noise is the Hurst exponent H. There is a loss of multifractality in combustion noise as combustors progress towards combustion instability, which is reflected in a decline of H 7,12 . We compute H 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 12 , 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 H drops to a very low value when the instability is present, the performance as an early warning signal is not convincing. Table IV shows the prediction accuracy using an optimal threshold for the test data set. When using 2 to 4 cycles for the computation of H, the prediction accuracy is even lower. The main reason for this is that H fluctuates strongly and therefore generates many false alarms. As an example, Fig. 8 displays the time development of H 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 safetycritical 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. autocorrelation function 51 and Cao's method 48 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 48 . 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.
By using recurrence plots analogous to the shown ones, RQA measures are calculated using the pyunicorn Python package 52 . 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 scikitlearn 53 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 C = 0.195 and a radial basis function kernel coefficient γ = 0.918.