1 Nonlinear time-delayed thermoacoustic model
We investigate the acoustics of a resonator excited by a heat source, which is a monopole source of sound. The main assumptions of the reduced-order model are : (i) the acoustics evolve on top of a uniform mean flow; (ii) the mean-flow Mach number is negligible, therefore the acoustics are linear and no flow inhomogeneities are convected; (iii) the flow is isentropic except at the heat-source location; (iv) the length of the duct is sufficiently larger than the diameter, such that the cut-on frequency is high, i.e., only longitudinal acoustics are considered; (v) the heat source is compact, i.e., it excites the acoustics at a specific location, ; (vi) the boundary conditions are ideal and open-ended, i.e., the acoustic pressure at the ends is zero. Under these assumptions, the non-dimensional momentum and energy equations read, respectively 
where is the acoustic velocity; is the acoustic pressure; is the time; is the axial coordinate of the duct; is the Dirac delta distribution at the heat source location, ; is the damping factor, which models the acoustic energy radiation from the boundaries and thermo-viscous losses; and is the heat release rate (or, simply, heat release). The heat release, , is modelled by a nonlinear time delayed law 
where is the time delay; is the strength of the heat source; and . The time delay and strength of the heat source are the two key parameters of a reduced-order model for the flame . Physically, is the time that the heat release takes to respond to a velocity perturbation at the flame’s base; while provides the strength of the coupling between the heat source and the acoustics. Velocity, pressure, length and time are nondimensionalized as in 
. The set of nonlinear time-delayed partial differential equations (1)-(2) provides a physics-based reduced-order model for the nonlinear thermoacoustic dynamics. Owing to the assumptions we made, the model can only qualitatively replicate the nonlinear thermoacoustic behaviour. In this paper, we propose a Lagrangian method to make a qualitative model quantitatively more accurate any time that reference data can be assimilated. Such data can come, for example, from sensors in experiments or time series from high-fidelity simulations.
1.1 Numerical discretization with acoustic modes
We use separation of variables to decouple the time and spatial dependencies of the solution. The spatial basis on to which the solution is projected consists of the natural acoustic modes. When decomposed on the natural acoustic eigenfunctions, the acoustic velocity and pressure read, respectively
where the relationship between and has yet to be found and is the number of acoustic modes considered. This discretization is sometimes known as the Galerkin discretization . The state of the system is provided by the amplitudes of the Galerkin modes that represent velocity, , and those that represent pressure, . The damping has modal components, , where and are damping coefficients [6, 9, 1, 5, 8]
. In vector notation,and . The state vector of the discretized system is the column vector . The governing equations of the Galerkin modes are found by substituting (4)-(5) into (1)-(3). Equation (2) is then multiplied by and integrated over the domain (projection). In so doing, the spatial dependency is removed and the Galerkin amplitudes are governed by a set of nonlinear time-delayed differential equations
where ; and . The labels are introduced for the definition of the Lagrangian (Sec. 3.2). Because the Galerkin expansions (4)-(5) are truncated at the -th mode, we obtain a reduced-order model of the original thermoacoustic system (1)-(2) with degrees of freedom (6)-(8). The reduced-order model is physically a set of time-delayed oscillators, which are nonlinearly coupled through the heat release term. In the following sections, we employ 4D-Var data assimilation to improve the accuracy of such a reduced-order model111Although different from our study, it is worth mentioning that a study that combined 4D-Var data assimilation with reduced-order models of the Navier-Stokes equations based on proper orthogonal decomposition can be found in ..
2 Data assimilation as a constrained optimization problem
The ingredients of data assimilation are (i) a reduced-order model to predict the amplitude of thermoacoustic oscillations, which provides the so-called background state vector at any time, (red thick line in Fig. 1); (ii) data from external observations, , which is available from high-fidelity simulations or experiments at times (grey diamonds in Fig. 1); and (iii) an educated guess on the parameters of the system, which originates from past experience. The uncertainties on the background solution and observations are here assumed normal and unbiased. and are the background and observation covariance matrices, respectively, which need to be prescribed. For simplicity, we assume that and
are diagonal with variancesand (i.e., errors are statistically independent). A cost functional is defined to measure the statistical distance between the background predictions and the evidence from observations. First, we want the state of the system to be as close as possible to the observations. Second, we do not want the improved solution to be “too far” away from the background solution. This is because we trust that our reduced-order model provides a reasonable solution. Mathematically, these two requirements can be met, respectively, by minimizing the following cost functional
The objective of state estimation is to improve the prediction of the state, x, over the interval , by reinitializing the background initial conditions, , with optimal initial conditions. These optimal initial conditions are called analysis initial conditions, , which are the minimizers of the cost functional (9). We start from a background knowledge of the model’s initial conditions, , which is provided by the reduced-order model when data is not assimilated. By integrating the system from , we obtain the red trajectory in Fig. 1, . The set of observations is assumed to be distributed over an assimilation window at some time instants. Pictorially, the analysis trajectory corresponds to the green thin line in Fig. 1, which is the minimal statistical distance between the background initial condition (magenta thick arrow) and observations (blue thin arrows). This algorithm is known as 4D-Var in weather forecasting . State estimation enables the adaptive update of reduced-order thermoacoustic models whenever data is available.
3 Data assimilation for nonlinear thermoacoustic dynamics
We propose a set of cost functionals to perform data assimilation with the thermoacoustic model introduced in Sec. 1. We also introduce the formalism to perform Lagrangian optimization, in which adjoint equations enable the efficient calculation of the gradients of the thermoacoustic cost functionals with respect to the initial state.
3.1 Thermoacoustic cost functionals
Crucial to data assimilation is the definition of the cost functionals and . Three physical thermoacoustic cost functionals are proposed and compared to reproduce different scenarios. For the background error
For the observation error
where is the location where the measurement is taken and is the instant at which the th observation is assimilated. On the one hand, by using and the analysis solution is optimized against the background pressure at and the measured pressure at , . Physically, this means that the acoustic pressure is the model output, , and the observable from the sensors, . On the other hand, and constrain the pressure modes. Physically, this means that every pressure mode provided by the background solution is a model output, , and it is assumed that the modes of the acoustic pressure,, can be calculated from the sensors on the fly. For the background cost functional, we also defined as a norm of the modes, which does not have a corresponding observation cost functional because the spatial dependency is not explicit.
To attain a minimum of , a necessary condition is that the gradient vanishes, i.e.,
where is the gradient with respect to the initial conditions. There exists such that because of the convexity of the cost functionals in the neighbourhood of a local extremum. To compute and , we use calculus of variation (Sec. 3.2). The Lagrange multipliers, also known as adjoint, or dual, or co-state variables (Sec. 3.3), provide the gradient information with respect to the initial state.
3.2 Lagrangian of the thermoacoustic system
The governing equations and their initial conditions are rewritten in the form of constraints, F, which hold over time intervals, while G are the constraints that hold for a specific time only, i.e., . Together with equations (6)-(8) and by defining the auxiliary variable , they read
The constraints for the initial conditions read
By defining an inner product
where and are arbitrary functions, the Lagrangian of the nonlinear system can be written as
where each is
where , , and are the Lagrange multipliers, or adjoint variables, of the corresponding constraints. Because we wish to derive the adjoint equations for the cost functional , we consider the time window to be .
3.3 Adjoint equations
We briefly report the steps to derive the evolution equations of the Lagrange multipliers (adjoint equations) . First, the Lagrangian (21) is integrated by parts to make the dependence on the direct variables explicit. Secondly, the first variation is calculated by a Fréchet derivative
Thirdly, the derivatives of (21) are taken with respect to the initial condition of each variable of the system, . These expressions will be used later to compute the gradient. Finally, to find the set of Lagrange multipliers that characterizes an extremum of the Lagrangian, , variations with respect to are set to zero. The adjoint equations and their initial conditions are derived by setting variations of , and to zero over .
3.4 Gradient-based optimization
The optimization loop consists of the following steps:
Integrate the system forward from to from an initial state ;
Initialize the adjoint variables;
Evolve the adjoint variables backward from to ;
Evaluate the gradient using the adjoint variables at .
Once the gradient is numerically computed, the cost functional can be minimized via a gradient based optimization loop. The conjugate gradient  is used to update the cost functional until the condition is attained to a relative numerical tolerance of . By using a gradient based approach, we find a local minimum of . We verify that there is no other local minimum by computing in the vicinity of .
We validate the data-assimilation algorithm by twin experiments: The true state solution, , is produced by perturbing the unstable fixed point at the origin of the phase space222Strong constraint 4D-Var assumes that the model is perfect and the uncertainty is only in the initial conditions, therefore, the true trajectory can be a model output.; the background trajectory, , is obtained by perturbing each true mode initial condition with Gaussian error with variance ; the th observation is produced by adding Gaussian error with variance to . The outcome of twin experiments is summarized by the error plots shown in Figs. 3 and 4. The cyan vertical line indicates the end of the assimilation window, the red thick line is the difference between the true pressure and the background pressure, the green thin line is the difference between the true pressure and the analysis pressure. First, it is shown how the number of computed acoustic modes affects the solution of the system. Secondly, we investigate the effects that the different cost functionals have on the analysis solution. Finally, we discuss the effects that the number of observations have on the analysis.
The parameters we use are , , , and for the heat release term . The position of the heat source is and all the measurements are taken at .
4.1 Remarks on the thermoacoustic nonlinear dynamics
The thermoacoustic model is a set of nonlinearly coupled oscillators, which we initialize by imposing non-equilibrium initial conditions. We compare two solutions, using and in Figs. 2a and 2b, respectively. Higher modes are quickly damped out, thus, after a transient where strong nonlinear modal coupling occurs, the solution obtained with is qualitatively similar to the solution obtained with . During the transient, if sufficient modes are computed, the dynamics are more unpredictable because of the intricate modal interaction. The twin experiments are performed with , which provide a more accurate solution. It is shown that state estimation is markedly affected depending on whether we observe the system during the transient or at regime.
4.1.1 Effect of the observation error
We can simulate two main scenarios, depending on the choice of (Sec. 3.1). Figure 3a is obtained using , i.e., by modelling observations about the pressure only. The analysis pressure error slightly deviates from zero in the assimilation window. When the forecast window starts, the analysis suddenly approaches the background again. Figure 3b is obtained using , therefore the observations contain information about every pressure modes. The forecast quality is considerably enhanced. The assimilation window is , thus, the observations are obtained during the transient, which lasts up to , where the dynamics are more unpredictable due to nonlinear interactions between modes. As we will show in the next subsection, increasing is not an effective strategy to improve the forecast during the transient when the pressure is observed.
4.1.2 Effect of the background error
From a numerical standpoint, the background error, , acts as an observation at . The analysis trajectory is an optimal blending of the information from the measurements and the previous educated model output. Therefore, we emphasize that the outcome of twin experiments is improved if the cost functionals of the background knowledge and observations are consistent. In the present framework, it means that should be used with and should be used with . On the one hand, if the source of assimilated data favours the background knowledge (e.g. using or together with ), the analysis trajectory will be closer to the background. On the other hand, if the source of assimilated data favours the observations (e.g. using with ), the analysis trajectory will be closer to the observations.
4.1.3 Effect of the number of observations
Generally speaking, the higher the number of observations the more the optimal solution will be similar to the true solution. This can be deduced by inspection of Figs. 4a and 4b. The value of is increased from 50 to 250, over an assimilation window of 2.5 time units (the observations are not shown in the figures), resulting in a smaller error amplitude when more observations are available.
However, it is possible that the increasing of will not result in a better state estimation. This happens if we measure the pressure during the transient, that is, using and . In the transient, the measured pressure is a combination of all modes. Therefore, the same pressure level is associated with different combinations of modes, hence, no useful information is assimilated to help determine the state of the system by knowing only the pressure (i.e., by using ). Under these circumstances, the forecast quality will remain poor, as shown in Fig. 3a, regardless of the number of observations. When the assimilation window is extended up to 2.5 time units, as shown in Fig. 4, the pressure is observed also at regime, that is in , approximately. In this interval, the measured pressure signal is produced chiefly by the first three modes because higher modes are damped out. Therefore, the measured pressure becomes more effective information about the state to assimilate. At regime, as intuitively expected, increasing the observation frequency produces a better forecast.
We conclude that having poor information about the system’s state cannot be simply balanced by increasing the number of observations. Given a number of observations with their time distrbution, it is the synergy between an appropriate cost functional and the recognition of the type of dynamics (transient vs. regime) to be the key for a successful data assimilation.
5 Conclusions and future work
Preliminary thermoacoustic design is based on simplified and computationally cheap reduced-order models that capture the inception of thermoacoustic instabilities (linear analysis) and their saturation to finite amplitude oscillations (nonlinear analysis). We propose a Lagrangian method to make a qualitative reduced-order model quantitatively more accurate any time that reference data can be assimilated. Such data can come, for example, from sensors in experiments or time series from high-fidelity simulations. To test the method we perform a series of twin experiments with the thermoacoustic model of a resonator excited by a heat source (horizontal Rijke tube). When sufficient modes are computed, a clear distinction emerges between a transient state and the state at regime. The former is characterized by irregular dynamics due to the interaction between all modes, while at regime the dynamics are chiefly dominated by the first three modes. We find that, at regime, it is possible to enhance the forecast by assimilating data about the pressure. As intuitively expected, the higher the number of observations, the better the forecast accuracy. While testing the effectiveness of data assimilation during the transient, we find that it is not possible to improve the forecast by measuring the pressure only. Moreover, the quality of the forecast remains poor regardless of the number of observations. Therefore, we propose a more effective cost functional, which takes into consideration the spectral content of the measured signal to enable a successful state estimation also during the transient. In state estimation, we implicitly assume that the parameters are correct. However, this is rarely the case in thermoacoustics, where the parameters are uncertain and need to be optimally calibrated. Ongoing work includes simultaneous parameter and state estimation using Lagrangian optimization with state augmentation. This work opens up new possibilities for on-the-fly optimal calibration and state estimation of reduced-order models in thermoacoustics for applications in propulsion and power generation.
-  Balasubramanian, K., Sujith, R.I.: Thermoacoustic instability in a Rijke tube: Non-normality and nonlinearity. Physics of Fluids 20(4), 044103 (2008). https://doi.org/10.1063/1.2895634
-  Blayo, É., Bocquet, M., Cosme, E., Cugliandolo, L.F.: Advanced Data Assimilation for Geosciences. Oxford University Press, first edn. (2015). https://doi.org/10.1093/acprof:oso/9780198723844.001.0001
-  Crocco, L.: Research on combustion instability in liquid propellant rockets. Symposium (International) on Combustion 12(1), 85–99 (1969). https://doi.org/10.1016/S0082-0784(69)80394-2
-  Du, J., Navon, I.M., Zhu, J., Fang, F., Alekseev, A.K.: Reduced order modeling based on POD of a parabolized Navier-Stokes equations model II: Trust region POD 4D VAR data assimilation. Computers and Mathematics with Applications 65(3), 380–394 (2013). https://doi.org/10.1016/j.camwa.2012.06.001
-  Juniper, M.P.: Triggering in the horizontal Rijke tube: non-normality, transient growth and bypass transition. Journal of Fluid Mechanics 667, 272–308 (nov 2011). https://doi.org/10.1017/S0022112010004453
-  Landau, L.D., Lifshitz, E.M.: Fluid Mechanics. Pergamon Press, second edn. (1987)
-  Magri, L.: Adjoint methods as design tools in thermoacoustics. Applied Mechanics Reviews 71(2), 020801 (2019). https://doi.org/10.1115/1.4042821
-  Magri, L., Juniper, M.P.: Sensitivity analysis of a time-delayed thermo-acoustic system via an adjoint-based approach. Journal of Fluid Mechanics 719, 183–202 (2013). https://doi.org/10.1017/jfm.2012.639
-  Matveev, K.I., Culick, F.E.C.: A model for combustion instability involving vortex shedding. Combustion Science and Technology 175, 1059–1083 (2003)
-  Orchini, A., Rigas, G., Juniper, M.P.: Weakly nonlinear analysis of thermoacoustic bifurcations in the Rijke tube. Journal of Fluid Mechanics 805, 523–550 (2016). https://doi.org/10.1017/jfm.2016.494
-  Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P.: Numerical recipes. Cambridge University Press, 3rd edn. (2007)
-  Zinn, B.T., Lores, M.E.: Application of the Galerkin Method in the Solution of Non-linear Axial Combustion Instability Problems in Liquid Rockets. Combustion Science and Technology 4(1), 269–278 (1971). https://doi.org/10.1080/00102207108952493