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 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 [20]. Those instabilities can be further subcategorized into highfrequency (screeching) and lowfrequency 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 highfrequency 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 timeseries 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 01 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 highamplitude 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 handcrafted 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 entropycomplexity 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
[43]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
[29]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 datadriven 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 stateoftheart 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 datadriven 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 selfexcited highfrequency 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 highfrequency combustion instabilities [27]. However, the current study only focuses on the prediction of thermoacoustic instabilities for the propellant combination LOX/H in BKD.
The BKD setup is shown in Fig. 1. It consists of three main components: a multielement injector head, a cylindrical combustion chamber, and a convergentdivergent 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 watercooled 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 watercooled highfrequency piezoelectric pressure sensors. Eight Kistler type 6043A are flushmounted in the ring with an even circumferential distribution, in order to measure the chamber pressure oscillations . The highfrequency pressure sensors have a measurement range set to bar and a sampling rate of 100 kHz. An antialiasing filter was set at 30 kHz.
Two different types of selfexcited highfrequency 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.
Stable and unstable operating conditions can be identified in the spectrogram. The strongest selfexcited highfrequency 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 peaktopeak amplitudes during unstable combustion reached up to 1635 bar but showed a large temporal variation. In addition to the pressure oscillation data, Gröning et al. [10, 12] used fibreoptical 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 injectiondriven [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 highspeed 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 LOXcoupled 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.
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 highamplitude 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 2630 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 threeway 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 LOXposts 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 nonnormal 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:
(1) 
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:(2) 
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 twodimensional 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:
(3) 
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 :
(4) 
where
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:(5) 
where
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 :(6) 
Using , one can further calculate the associated Shannon entropy (ENTR):
(7) 
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(8) 
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
(9) 
where and . The associated function is of the form
(10) 
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:
(11) 
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 subgradient 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,
(12) 
and the decision function can be written as
(13) 
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
(14) 
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
(15) 
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,
(16) 
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 (limitcycle) 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 datadriven 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 state^{1}^{1}1GH (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 highpressure 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 length^{2}^{2}2Most experiments used a LOXpost 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 
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 peaktopeak 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 peaktopeak 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 peaktopeak 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 timedependent 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 crossvalidation to assess the performance of the SVM measured by the Fscore (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 leadtime 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 truepositive 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 falsepositive 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.
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 % 
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.
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 Fscore. 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 % 
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 % 
Viii Conclusion and Outlook
In this paper, we presented the datadriven construction of a warning signal for the early detection of thermoacoustic instabilities in rocket thrust chambers. Using timeseries data from experimental tests performed with the cryogenic rocket thrust chamber BKD, we calculated timedependent 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, highfrequency 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, socalled outofdistribution 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 resourceintensive 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.
Acknowledgements.
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 highpass filter is used to remove nonacoustic 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 timedependent 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 [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 scikitlearn
[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 constantand a radial basis function kernel coefficient
.Appendix B BKD Runs
References
 [1] (2018) Experimental investigation of selfexcited 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] (2019) Injectordriven flame dynamics in a highpressure multielement oxygenhydrogen rocket thrust chamber. Journal of Propulsion and Power 35, pp. 632–644. External Links: Document Cited by: §II, §II.
 [3] (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] (2006) Pattern recognition and machine learning. Springer. Cited by: §IV, §V.
 [5] (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] (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] (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] (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] (2011) Dynamic properties of combustion instability in a lean premixed gasturbine combustor. Chaos. External Links: Document, ISSN 10541500 Cited by: §I.
 [10] (2016) Injectordriven 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] (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] (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] (2017) Untersuchung selbsterregter verbrennungsinstabilitäten in einer raketenbrennkammer. Ph.D. Thesis, RWTH Aachen. Cited by: §II, §II, §II.
 [14] (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] (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] (2015) Large eddy simulations of multiple transcritical coaxial flames submitted to a highfrequency transverse acoustic modulation. Proceedings of the Combustion Institute 35 (2), pp. 1461–1468. External Links: Document Cited by: §II.
 [17] (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] (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 SP194. Cited by: §I, §I.
 [20] (2019) Rocket propulsion. 1 edition, Cambridge University Press. Cited by: §I.
 [21] (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] (2018) Sensitivity and nonlinearity of thermoacoustic oscillations. Annual Review of Fluid Mechanics. External Links: Document Cited by: §I, §III, §V, §VII.
 [23] (2011) Triggering in the horizontal rijke tube: nonnormality, transient growth and bypass transition. Journal of Fluid Mechanics. External Links: Document Cited by: §II.
 [24] (2019) Dynamical systems approach to study thermoacoustic transitions in a liquid rocket combustor. Chaos. External Links: Document Cited by: §III, §V.
 [25] (201906) 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] (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] (2020) Experimental investigation of selfexcited combustion instabilities in a lox/lng rocket combustor. In AIAA Scitech 2020 Forum, AIAA Paper 20200423, Orlando, Florida. External Links: Document Cited by: §II.
 [28] (2003) Encounters with neighbours: current developments of concepts based on recurrence plots and their applications. Norbert Marwan. Cited by: §III.
 [29] (2020) Online prediction of combustion instabilities using machine learning. Journal of Engineering for Gas Turbines and Power. Cited by: §I.
 [30] (2018) Early detection of thermoacoustic instabilities using hidden markov models. Combustion Science and Technology. Cited by: §I.
 [31] (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] (201405) Multifractality in combustion noise: predicting an impending combustion instability. Journal of Fluid Mechanics 747, pp. 635–655. External Links: Document, ISSN 00221120, Link Cited by: §I, §VII.
 [33] (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/17568277.5.4.273 Cited by: §I.
 [34] (2014) Intermittency route to thermoacoustic instability in turbulent combustors. Journal of Fluid Mechanics 756, pp. 470. Cited by: §I.
 [35] (2008) Applied nonlinear dynamics: analytical, computational, and experimental methods. John Wiley & Sons. Cited by: Appendix A.
 [36] (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 20116027, San Diego, California. External Links: Document Cited by: §II.
 [37] (2011) Scikitlearn: machine learning in Python. Journal of Machine Learning Research 12, pp. 2825–2830. Cited by: Appendix A.
 [38] (1878) The explanation of certain acoustical phenomena. Nature 18 (455), pp. 319–321. Cited by: §I.
 [39] (2016) Dynamic datadriven prediction of instability in a swirlstabilized combustor. International Journal of Spray and Combustion Dynamics 8 (4), pp. 235–253. Cited by: §I.
 [40] (2017) Largeeddy simulations of a subscale liquid rocket combustor: influence of fuel injection temperature on thermoacoustic stability. In 7th European Conference for Aeronautics and Aerospace Sciences (EUCASS) 2017, Milano, Italy. External Links: Document Cited by: §II.
 [41] (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] (2016) "L42" technology demonstrator: operational experience. In Space Propulsion 2016, Rome, Italy. Cited by: §II.
 [43] (2020) Bayesian machine learning for the prognosis of combustion instabilities from noise. Journal of Engineering for Gas Turbines and Power. Cited by: §I.
 [44] (1980) Detecting strange attractors in turbulence. In Dynamical Systems and Turbulence, Warwick 1980, Warwick, United Kingdom. Cited by: §III.
 [45] (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]
(1999)
An overview of statistical learning theory
. IEEE Transactions on Neural Networks 10 (5), pp. 988–999. External Links: Document Cited by: §IV.  [47] (1966) Effect of propellant injection velocity on screech in 20000pound hydrogenoxygen rocket engine. Technical report NASA TN D3373, National Aeronautics and Space Administration, National Aeronautics and Space Administration, Cleveland, Ohio (USA). Lewis Research Center, Washington D.C.. Cited by: §II.
 [48] (202004) 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 15336808, Document Cited by: §I.

[49]
(2020)
A Reinforcement Learning Approach for Transient Control of Liquid Rocket Engines
. arXiv:2006.11108 [cs, eess, math, stat]. Cited by: §I.  [50] (2020) HardwareInTheLoop Tests of Complex Control Software for Rocket Propulsion Systems. In Proceedings of the 71st International Astronautical Congress, Virtual. Cited by: §I.
 [51] (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] (2005) Combustion instabilities: basic concepts. Combustion Instabilities in Gas Turbine Engines. Cited by: §II.
Comments
There are no comments yet.