Neural ODE to model and prognose thermoacoustic instability

Figure 1 shows the feed forward NN architecture used in this study. NN maps the input data to the output through series of mathematical manipulations. The layer of nodes that receives the input is the input layer and is denoted as Layer ($0$). There are $n^{(0)}$ inputs represented as vector ${Z}^{(0)}=[{z}^{(0)}_{1},..,{z}^{(0)}_{n^{(0)}}]^{T}$. The layer of nodes that outputs the desired objective is Layer ($k$). There are $n^{(k)}$ outputs represented by the vector, ${Z}^{(k)}=[{z}^{(k)}_{1},..,{z}^{(k)}_{n^{(k)}}]^{T}$. The layers of nodes that lie between Layer ($0$) and Layer ($k$); i.e., layers $(1)$ to $(k-1)$, are called hidden layers. Biological neurons respond or generates action potential if the aggregate input stimuli is above a certain threshold. A similar behavior can be approximated in artificial NN at each node. A node maps the weighted sum of inputs to its output using a nonlinear mapping function called the activation function ($A$). There are different types of activation functions such as $sigmoid$, $tanh$, $ReLu$ Szandała (2021). Thus the output of the node $(j)$ in the Layer $(l)$; i.e., $z^{(l)}_{j}$ is computed as,

where, $w^{(l)}_{ij}$ is the weight assigned to the link that connects node $(i)$ in Layer $(l-1)$ to node $(j)$ in layer $(l)$ and $n^{(l-1)}$ is the number of nodes in Layer $(l-1)$. The vector of outputs of all nodes in layer $(l)$ is, $Z^{(l)}=[z^{(l)}_{1},..,z^{(l)}_{n^{(l)}}]^{T}$. Thus, $Z^{(l)}$ can be represented as,

where, $W^{(l)}=(w_{ij})\in\mathbb{R}^{n_{(l-1)}\times{n_{(l)}}}$, is the weight matrix. Defining a NN involves specifying the number of layers, nodes in each layer, activation function and weight matrices. This modular structure allows to modify the complexity of the function described by NN.

We denote the augmented matrix of all the weight matrices as, $\theta=[W^{(1)},..,W^{(k)}]$. Thus, $\theta$ contains all the parameters of the NN. The NN can be expressed as a function ($f$) that maps the input vector $Z^{(0)}$ to the output vector $Z^{(k)}$,

For further discussion, we treat $X=Z^{(0)}$ as the input and $Y=Z^{(k)}$ as the output predicted by the NN. Thus,

The data set obtained from experiment is, $\mathcal{D}=\{(X_{i},Y_{tr,i})|i=(1,..,N)\}$, consisting of $N$ data samples, where $Y_{tr,i}$ is the true measured output from the experiment for the input $X_{i}$. Now, consider a NN which produces $Y_{i}$ as the output when the input is $X_{i}$. The penalty for failing to predict true $Y_{tr,i}$ and instead predicting $Y_{i}$ can be quantified using a loss function $L$,

The average loss over the complete data set is,

Training a NN involves optimizing the loss function $\mathcal{L}$ with respect to $\theta$. Based on the training objective, the appropriate loss function can be constructed Janocha and Czarnecki (2016); Bishop (2006). The weight matrices corresponding to the different layers of the NN are randomly initialized before training. The process of optimization can be performed by gradient descend which involves marching the weights in the opposite direction of the gradient of the loss function with respect to $\theta$ Curry (1944). The gradient is often computed using the back-propagation algorithm Rumelhart, Hinton, and Williams (1986). Considering $\theta_{n}$ as the value of $\theta$ at $(n)^{th}$ iteration, during the $(n+1)^{th}$ iteration, the gradient decent algorithm updates this value of $\theta$ to $\theta_{n+1}$ as,

where, $\alpha$ is the learning rate. Training a NN involves minimizing the loss function with iterative application of this update rule. The discrepancy between the true (expected) output $Y_{tr}$ and the predicted output $Y_{pred}$ reduces upon repeated application of the update rule. Thus, trained or optimized NN can approximately predict the output.

A dynamical system that is continuous in time can be represented in the form of a differential equation. NNs can be trained to approximate the differential equation of the dynamical system generating the data. This can be achieved by training NN to output the time derivative of the input state vector. As discussed before, the input vector as $X$ to the NN and the output vector as $Y$. This means that the output of NN, representing differential equation of the dynamical system, should output $Y$ which is approximately equal to time derivative of input state vector ${X}$; i.e., $Y\approx\dot{X}$. A NN that approximates the ordinary differential equation (ODE) of the dynamical system is termed as a neural ODE Chen et al. (2018). Hence, the functional form of a neural ODE is,

The architecture of the neural ODE is the same as the previously discussed NN. Due to the constraint that the output is the time derivative of the input, the number of nodes in the output layer is the same as that of the input layer nodes. In Fig.1, variables in the braces of output layer forms the output vector of the neural ODE; i.e. $Y=[\dot{z}_{1}^{(0)},\dot{z}_{2}^{(0)},...,\dot{z}_{n^{(k)}}^{(0)}]=\dot{X}$. The evolution of the dynamical system with time can be obtained by integrating the neural ODE numerically. When $X_{tr}(t_{0})$ is the true initial state of the system at time $t_{0}$ (known from the data), the neural ODE predicts the the state of the system as $X(t_{i})$ at time $t_{i}$ as,

If, $X_{tr}(t_{i})$ is the true state of the system at time $t_{i}$, then the loss incurred due to prediction $X(t_{i})$ is,

Optimization of the loss function with respect to $\theta$ can be performed using the adjoint optimization method, details of which can be found in Chen et al. (2018). The parameter vector $\theta$ obtained after training is characteristic to the dynamical state of the system.

A set of variables $(x_{1},x_{2},..,x_{d})$ that gives the complete description of the state of system are known as state variables. Here, $x_{i}$ is a state variable and there are total $d$ state variables. The $d$ dimensional space with state variables as coordinates is referred to as the phase space. Thus, the phase space contains all possible states that the system can assume. In practical systems only a finite number of state variables can be measured. Under this condition, in order to visualize and analyze the asymptotic dynamics of an experimental system, we can reconstruct a phase space according to delay embedding theorem proposed by Takens (1981). For thermoacoustic systems the unsteady acoustic pressure, $p^{\prime}(t)$ can easily be measured. The evolution of the thermoacoustic system can be visualized by reconstructing the phase space using delay-coordinates. The delay coordinate can be expressed as, $[p^{\prime}(t),p^{\prime}(t+\tau),..,p^{\prime}(t+(d-1)\tau)]$. Here, $\tau$ is the delay and $d$ is the dimension of the reconstructed phase space known as the embedding dimension. The dimension of the reconstructed phase space should be sufficiently high to avoid overlap of the trajectories and false neighbors. Reconstruction of phase space involves two steps 1) finding the optimal time delay $\tau$ and 2) finding the dimension of the phase space Abarbanel et al. (1993); Kantz and Schreiber (2003). For the present analysis we used the method of average mutual information (AMI) to obtain $\tau$. The method has its roots in the information theory. The optimal delay is chosen such that the original time series $[p^{\prime}(t_{1}),p^{\prime}(t_{2}),..]$ conveys the least information about the delay time series $[p^{\prime}(t_{1}+\tau),p^{\prime}(t_{2}+\tau),..]$.

We used the method of false nearest neighbors (FNN) to obtain the embedding dimension. The co-ordinates of the dynamical system at a time instance $t_{i}$ with embedding dimension $d$ is $y_{d}(t_{i})=[p^{\prime}(t_{i}),p^{\prime}(t_{i}+\tau),..,p^{\prime}(t_{i}+(d-1)\tau)]$. The $d$ for which the total number of FNN drops to zero is chosen as the embedding dimension. A brief description of AMI and FNN is given in Appendix A.1. The state vector for the reconstructed phase space using the time series of pressure is, $X(t)=[p^{\prime}(t),p^{\prime}(t+\tau),..,p^{\prime}(t+(d-1)\tau)]$. The dynamical system is then represented as:

This equation is the same as equation (8), where ${f}$ is the function representing a NN. We construct NN with $tanh$ activation function and one hidden layer. According to universal approximation theorem, a NN even with a single hidden layer with a smooth activation function can approximate any function if sufficiently hidden units or nodes are considered Kurt, Maxwell, and White (1989). The $tanh$ activation function is used as it is smooth, changes sign across the zero and could approximate smooth functions better compared to piecewise linear functions Szandała (2021). The absolute value of the error is used to construct the loss function ($L$) as:

During training, the absolute value of the average loss function is optimized with respect to the neural ODE parameter vector $\theta$ at a fixed Reynolds number ($Re$). Thus the parameters of the neural ODE depends on the Reynolds number; i.e., $\theta=\theta(Re)$. The neural ODE trained on the experimental data gives an approximate dynamical system representation for the experimental system.

Thermoacoustic systems are characterized by the mutual interaction of the heat source and the acoustic field in the confinement. Pressure fluctuations inside the combustor can easily be monitored using a microphone. Measurements of heat release rate fluctuations are often unavailable due to inaccessibility of the combustor to photo-multiplier tube. Thus, one should be able to predict the onset of periodic oscillations using the time series of pressure fluctuations alone. We reconstruct the phase space using the time series of pressure as, $X_{p}(t)=[p^{\prime}(t),p^{\prime}(t+\tau_{p}),..,p^{\prime}(t+(d_{p}-1)\tau_{p})]$. Here, $\tau_{p}$ and $d_{p}$ are the time delay and embedding dimension respectively, for the time series of pressure fluctuations. The dynamical system is then represented as:

The construction of loss function and its optimization are already discussed in II.3. The present model constructed based on the time series of pressure fluctuations has parameter vector $\theta$ that governs the behavior of the model. A smooth variation of these parameters can lead to a qualitative change in the dynamical state of the system known as bifurcation.

The time series of $p^{\prime}$ with $N$ data points is measured for each operating state defined by the $Re$ of the flow. We train our NN on the first $n_{w}$ continuous data samples. After training with this window of $n_{w}$ data points, the window is shifted by $n_{s}$ data points and the neural ODE is trained again with the new $n_{w}$ data points. Therefore, for each operating state, $[N/n_{s}]$ different neural ODEs and their corresponding parameter vectors $\theta$ are obtained ([.] is the greatest integer function). To get single representative $\theta$ for a given $Re$, we average the parameter vector to get $\bar{\theta}(Re)=[\overline{{W}}^{(1)},..,\overline{{W}}^{(k)}]$.

A sine wave with a frequency ($f$ Hz) equal to the limit cycle oscillations ($f_{LC}$ Hz) is constructed (we observed that for 5% deviation in $f$ from $f_{LC}$, the proposed precursor measures shows the expected behavior). The neural ODE is trained on this sine wave and the obtained parameter vector, $\theta_{ref}=[{W}^{(1)}_{ref},..,{W}^{(k)}_{ref}]$ is used as a reference. The cosine of the angle between the ${{W}}^{(l)}(Re)$ and ${{W}}^{(l)}_{ref}$ quantifies the alignment of the weight matrix of the current operating state at specified $Re$ with the reference weight matrix obtained from sine wave. This precursor measure $m(Re)$ is obtained as,

where, $<,>$ denotes the inner product and $|.|_{2}$ denotes the $L_{2}$ norm of the matrix.

For turbulent combustors the reacting flow field is the source of heat. Modeling a thermoacoustic system using linearized wave equation needs a model for the heat source. The response of the heat release rate to velocity fluctuations is often modeled using a flame transfer function or a flame describing function. The flame transfer function models the heat release rate as a linear response to the velocity fluctuations and hence the gain and phase of the transfer function are independent of the amplitude of the acoustic velocity fluctuations. However, the flame response is a nonlinear function of the amplitude of the acoustic velocity fluctuations. To capture this nonlinear dependence, researchers use a flame transfer function with dependence on acoustic fluctuation amplitude known as the flame describing function. In general, the heat release rate is modulated by the oscillations in the flame surface, equivalence ratio, acoustic velocity and pressure fluctuations. The net heat release rate can be obtained by integrating the heat release rate over the flame surface.The heat release rate can, in general, expressed as,

where, ${u}^{\prime}$ and ${p}^{\prime}$ are the velocity and pressure fluctuations. In present study, we are using the net heat release rate ($\dot{q}$) and not the mean subtracted ($\dot{q}^{\prime}$).

In present study, the time series of pressure ($p^{\prime}$) and heat release rate fluctuations ($\dot{q}$) are obtained from experiments conducted on a turbulent combustor. Measurements of the velocity fluctuations $u^{\prime}$ are absent. As discussed before, the reconstruction of phase space can be used to create a high dimensional phase space to describe the dynamics. Thus, dynamical system representation can be obtained for the thermoacoustic system, even in the absence of measurement of $u^{\prime}$. To reconstruct the phase space, the delay and the embedding dimension for the $p^{\prime}$ and $\dot{q}$ are individually obtained as ($\tau_{p}$, $d_{p}$) and ($\tau_{q}$, $d_{q}$). The state vector for the reconstruction of phase is constructed as, $X_{pq}$ = $[p^{\prime}(t),p^{\prime}(t+\tau_{p}),..,p^{\prime}(t+(d_{p}-1)\tau_{p}),\dot{q}(t),\dot{q}(t+\tau_{q}),..,\dot{q}(t+(d_{q}-1)\tau_{q})]$. Thus, thermoacoustic system can be represented as,

We use neural ODE to approximate this mapping of state vector to its time derivative. On training the neural ODE on time series of $p^{\prime}$ and $\dot{q}$, we get a model for the thermoacoustic system. The absolute value of the difference between the predicted and the true time series is used as loss function. Optimization of the average loss function with $\theta$ gives the parameter vector for the neural ODE. Since the neural ODE is trained on time series representing the simultaneous variation of $p^{\prime}$ and $\dot{q}$, a NN with the optimum $\theta$ can represent the coupled behavior of the unsteady heat release rate and acoustic perturbations in a thermoacoustic system.

A thermoacoustic system involves the interaction between the acoustic field and the heat the release rate. Thus, we use $X_{pq}$ as the state vector as discussed in Sec. II.5.

The NN considered for the present analysis, contains only a single hidden layer with $tanh$ activation function. For the present analysis, the hidden layer always contains 100 nodes (the convergence of the predicted time series to the experimental data is given in Appendix A.2). The input and output layer contains nodes equal to the dimension of the state vector.

To visualize the change in the dynamical state of the system, all these time series, obtained for a sequence of $Re$, are stacked and are shown in Fig. 3(a). We train the neural ODE on the first $n_{w}$ data points. To compare this length of data with the number of data points $n_{c}$ in one period of limit cycle oscillation, we use $N_{c}=n_{w}/n_{c}$. We used $N_{c}\approx 48$ to train the network. We shift the training window by $n_{s}$ data points to train another NN. Thus, for each $Re$ we have $[N/n_{s}]$ parameter vectors. As discussed, to get the single representative neural ODE parameter vector, the average parameter vector $\bar{\theta}(Re)$ is calculated. The training is performed with a fixed number of $2000$ iterations and a learning rate of $\alpha=10^{-2}$.

The phase space reconstruction is performed using $AMI$ and $FNN$ methods as discussed in Sec. II.5. The time series data of periodic oscillations is used to obtain $AMI$ and $FNN$. The variation of the $AMI$ with delay ($\tau$) and $FNN$ with embedding dimension $d$ are shown in Fig. 7 in Appendix A.1 . For the time series of $\dot{q}$ and $p^{\prime}$, the $AMI$ attains a first minima at $\tau_{q}$ = 22 and $\tau_{p}$ = 21 respectively; these values ($\tau_{q}$ and $\tau_{p}$) are chosen as optimal delay. For time series of $\dot{q}$ and $p^{\prime}$ the number of $FNN$ become zero at dimension $d_{q}=10$ and $d_{p}=7$ respectively. The combined dimension $d=d_{q}+d_{p}$ is chosen as the embedding dimension to reconstruct the phase space.

The phase space trajectory of thermoacoustic system is reconstructed in ($p^{\prime}(t)$, ${p}^{\prime}(t+\tau)$, ${p}^{\prime}(t+2\tau)$) coordinates for visualization. The phase space trajectory constructed from experimental data for the state of chaotic, intermittent and limit cycle oscillations are shown in Fig. 3 (e), (f) and (g) respectively. These time series are reconstructed with the neural ODE and are shown in Fig. 3 (h), (i) and (j) respectively. This reconstruction is performed by using the same parameter ${\theta}(Re)$ and the initial condition as the corresponding time series obtained from experiments shown in Fig. 3 (e), (f) and (g) respectively. Integration is performed with an explicit $4^{th}$ order Runge-Kutta integrator Ernst Hairer (1993) for the duration of 0.4 s. Considering time period of limit cycle oscillations as $T_{LC}$, we integrate the neural ODE for close to $48T_{LC}$ duration.

For the experimental data, the phase space trajectory for the chaotic state shows that there is a cluttered trajectory, occupying only a small region in the phase space as shown in Fig. 3 (e). This figure also shows that the trajectory does not settle to a fixed point. The same behavior can be observed for the reconstructed phase trajectory using the neural ODE shown in Fig. 3 (h). For the state of intermittency, Fig. 3 (i),the reconstructed phase space trajectories using neural ODE shows different orbits indicating switching of the system dynamics between chaos and periodic oscillations which is not clearly visible in Fig. 3 (f) due to overlapping of multiple bursts of periodic oscillations with slightly different amplitude. Fig. 3 (g) shows that for the state of limit cycle oscillations the trajectory forms a closed loop for experimental data and matches well with the reconstructed phase space trajectory obtained using neural ODE which is shown in Fig. 3 (j). Considering that the integration is performed for $48T_{LC}$ the phase space trajectories show that the neural ODE could predict the asymptotic behavior for the states of chaos, intermittency and limit cycle oscillations. Thus, Fig. 3 clearly shows that the neural ODE performs well on the training data.

The neural ODE is tested with initial condition outside the training data, and the comparison of time series is shown in Fig. 4. The neural ODE is integrated for a duration $\approx 12T_{LC}$. Figure 4 (a-c) and 4 (d-f) respectively compares the true and predicted time series of pressure $p^{\prime}$ and heat release rate $\dot{q}$ fluctuations for the state of chaos, intermittency and limit cycle oscillations. These dynamical states corresponds to the $Re$ = $1.07\times 10^{4}$, $1.32\times 10^{4}$, $1.92\times 10^{4}$ respectively. For chaos, the predicted trajectories for $p^{\prime}$ and $\dot{q}$ follows the true trajectory for $\approx 3T_{LC}\approx 3.3T_{\lambda}$ (here $T_{\lambda}$ is the Lyapunov time scale) which is marked by dashed vertical line in the figure and after that deviation grows. For the state of intermittency, the predicted and true trajectories for $p^{\prime}$ matches well for $\approx 8T_{LC}$ and indicated by dashed vertical line in the figure. After that the neural ODE qualitatively captures the true time series. For the state of limit cycle oscillations, comparison of predicted and true time series of $p^{\prime}$ and $\dot{q}$ shows that predictions agree well with true data. As the neural ODE predicts the time derivative of the state variables, Fig. 4 (g-i) and (j-l) shows the comparison of true and predicted time derivative of $p^{\prime}$ and $\dot{q}$ (i.e., $\dot{p}^{\prime}$ and $\ddot{q}$) for the state of chaos, intermittency and limit cycle oscillations, respectively. The time derivative of the experimental data is obtained using a $2^{nd}$ order central difference scheme. Ideally, if the true and predicted values match; i.e., $((.)_{true}=(.)_{pred})$ then the point should fall on the reference diagonal line as shown in Fig. 4 (g-l). For all the three states considered here, all the points fall in the close vicinity of the diagonal line. Thus, neural ODE captures the dynamics of the thermoacoustic system for all the three states and performs well over testing data.

For turbulent combustor, thermoacoustic instability is presaged by the formation of coherent structures, intermittency and a reduction in the dimensionality as discussed in Sec. I. For dynamical systems, smooth variation of the parameters of the differential equation can alter the dynamical state of the system. Such a bifurcation could also occur in the neural ODE on changing the parameter vector $\theta$. The dynamics of thermoacoustic system changes with $Re$; hence $\theta$ should be a function of $Re$ ; i.e., $\theta=\theta(Re)$. Now, we create a parametric space where the parameter vector $\theta$ is the position vector for the state of the thermoacoustic system. The position of the thermoacoustic system in parametric space changes on changing the $Re$. We choose ten elements from each weight matrix of the NN to visualize the evolution of elements of weight matrix with $Re$. The variation of these elements with $Re$ for $W^{(1)}$ is shown in the Fig. 5 (a) and for $W^{(2)}$ is shown in Fig. 5 (b). For the state of chaos with small $Re$, some elements of weight matrix $W^{(1)}$ shows high magnitude. As the thermoacoustic system approaches state of limit cycle oscillations by increasing $Re$, the magnitude of the weights approaches to zero. On the other hand weights from $W^{(2)}$ are close to zero for the state of chaos and the magnitude of the weights increases as the thermoacoustic system approaches limit cycle oscillations. This observed pattern shows that weight matrices hold the information about the proximity of the thermoacoustic system to the state of limit cycle oscillations.

We analyze the evolution of the thermoacoustic system in the parametric space with the elements of the parameter vector $\theta$ as coordinate. The point that describes the position of the thermoacoustic system in parametric space moves closer to the point that describes the reference state, obtained using sinusoidal signal, on increasing the $Re$. As the point describing the dynamic state of the system in the parametric space approaches the reference state, the parameter vector of the system becomes aligned to the reference vector. To quantify the proximity of the thermoacoustic system to the onset of limit cycle oscillations we propose a measure $m^{(l)}$ as discussed in Sec. II.4. As discussed before, we have $[N/n_{s}]$ parameter vectors for each operating state. On evaluating $m^{(l)}$ for all $[N/n_{s}]$ parameter vectors we compute the average measure $\bar{m}^{(l)}(Re)$ and the standard deviation $\sigma(Re)$. Thus, for the sequence of Reynolds numbers ($Re$) we have $\bar{m}^{(l)}(Re)$ and $\sigma(Re)$. The variation of $m^{(l)}$ with $Re$ is shown in Fig. 6. The marker shows the mean value of the measure $\bar{m}^{(l)}(Re)$ and $\pm\sigma(Re)$ bound is shown using the color band. The figure includes plots for variation of $m^{(1)}$ and $m^{(2)}$ that corresponds to the weight matrix $W^{(1)}$, $W^{(2)}$ respectively. Both the figures show that $m$ approaches a value of one as the system approaches the state of limit cycle oscillations. Thus, we propose that $m$ could be used as a measure to quantify the proximity of the dynamical state of the thermoacoustic system to the onset of thermoacoustic instability. To assess the generality of the proposed framework to detect the proximity of the system dynamics to oscillatory instability, we apply this methodology to an aeroacoustic system with turbulent flow as described in Appendix A.3.