Combined State and Parameter Estimation in Level-Set Methods

TL;DR

This paper introduces a Bayesian data assimilation framework using ensemble Kalman filters for level-set methods, enabling accurate state and parameter estimation in complex interface dynamics.

Contribution

It develops a novel physics-informed, data-driven Bayesian approach for joint state and parameter estimation in level-set models, validated through multiple test cases.

Findings

Full recovery of reference solutions in twin experiments

Effective uncertainty quantification in complex interface dynamics

Successful prediction of nonlinear flame behavior with topological changes

Abstract

Reduced-order models based on level-set methods are widely used tools to qualitatively capture and track the nonlinear dynamics of an interface. The aim of this paper is to develop a physics-informed, data-driven, statistically rigorous learning algorithm for state and parameter estimation with level-set methods. A Bayesian approach based on data assimilation is introduced. Data assimilation is enabled by the ensemble Kalman filter and smoother, which are used in their probabilistic formulations. The level-set data assimilation framework is verified in onedimensional and two-dimensional test cases, where state estimation, parameter estimation and uncertainty quantification are performed. The statistical performance of the proposed ensemble Kalman filter and smoother is quantified by twin experiments. In the twin experiments, the combined state and parameter estimation fully recovers the reference solution, which validates the proposed algorithm. The level-set data assimilation framework is then applied to the prediction of the nonlinear dynamics of a forced premixed flame, which exhibits the formation of sharp cusps and intricate topological changes, such as pinch-off events. The proposed physics-informed statistical learning algorithm opens up new possibilities for making reduced-order models of interfaces quantitatively predictive, any time that reference data is available.

Figures

Tables

Table 1. Table 1 : Conditional probability distributions in data assimilation.

Task	PDF	Description
State estimation	$p\left(x_{0:N}\mid y_{1:N},\theta ,f\right)$	Given a physical model and its parameters, what is our belief in a series of states?
Parameter estimation	$p\left(\theta \mid y_{1:N},f\right)$	Given a physical model, what is our belief in a set of parameters?
Model comparison	$p\left(f\mid y_{1:N}\right)$	Between two physical models, in which one do we believe more?

Table 2. Table 2 : Conditional probability distributions in state estimation.

Task	PDF	Description
Filtering	$p\left(x_k\mid y_{1:k},\theta ,f\right)$	Given all observations from the past and now, what is our belief in the current state?
Smoothing	$p\left(x_k\mid y_{1:N},\theta ,f\right)$	Given all observations from the past, the future and now, what is our belief in the current state?
Prediction	$p\left(x_{N+1}\mid y_{1:N},\theta ,f\right)$	Given all observations from the past and now, what is our belief in a future state?

Table 3. Table 3 : Conditional probability distributions in combined state and parameter estimation.

Task	PDF	Description
Filtering	$p(x_k,{\theta }_k\mid y_{1:k},f)$	Given all observations from the past and now, what is our belief in the current state and set of parameters?
Smoothing	$p(x_k,{\theta }_k\mid y_{1:N},f)$	Given all observations from the past, the future and now, what is our belief in the current state and set of parameters?

Table 4. Table 4 : G 𝐺 G -equation parameters for the reference simulation.

Parameter	$r_i$ [m]	$r_a$ [m]	$U$ [m/s]	$s_L$ [m/s]	$\omega$ [rad/s]	$K$ [-]	$\epsilon$ [-]
Parameter value	0.006	0.012	1.0	0.164	277.78	1.0	0.1

Equations

$$x_k=x_{k-1}+{\int }_{t_{k-1}}^{t_k}f(x(t),\theta )\mathrm{d}t=G(x_{k-1},\theta )\sim p(x_k\mid x_{k-1},\theta ,f),$$

$$y_k=M(x_k)\sim p(y_k\mid x_k).$$

$$p(x_k\mid x_{0:k-1},y_{1:k-1},\theta ,f)=p(x_k\mid x_{k-1},\theta ,f),$$

$$p(x_k\mid x_{k+1:N},y_{k+1:N},\theta ,f)=p(x_k\mid x_{k+1},\theta ,f),$$

$$p(y_k\mid x_{0:N},y_{1:k-1},y_{k+1:N},\theta ,f)=p(y_k\mid x_k).$$

$${\tilde{x}}_k=\left(\begin{array}{c}{x}_k\\ {\theta }_k\end{array}\right),\tilde{f}(\tilde{x}(t))=\left(\begin{array}{c}f(x(t),{\theta }_k)\\ 0\end{array}\right),{\tilde{x}}_k={\tilde{x}}_{k-1}+{\int }_{t_{k-1}}^{t_k}\tilde{f}(\tilde{x}(t))\mathrm{d}t.$$

$$p(x_k\mid y_{1:k},\theta ,f)$$

$$=\frac{p(y_k\mid x_k)p(x_k\mid y_{1:k-1},\theta ,f)}{p(y_k\mid y_{1:k-1},\theta ,f)}.$$

$$p(x_k\mid y_{1:k-1},\theta ,f)$$

$$=\int p(x_k\mid x_{k-1},y_{1:k-1},\theta ,f)p(x_{k-1}\mid y_{1:k-1},\theta ,f)\mathrm{d}{x}_{k-1}$$

$$=\int p(x_k\mid x_{k-1},\theta ,f)p(x_{k-1}\mid y_{1:k-1},\theta ,f)\mathrm{d}{x}_{k-1}.$$

$$p(x_k\mid y_{1:k-1},\theta ,f)=\int p(x_k\mid x_{k-1},\theta ,f)p(x_{k-1}\mid y_{1:k-1},\theta ,f)\mathrm{d}{x}_{k-1}.$$

$$p(x_k\mid y_{1:k},\theta ,f)=\frac{p(y_k\mid x_k)p(x_k\mid y_{1:k-1},\theta ,f)}{p(y_k\mid y_{1:k-1},\theta ,f)}.$$

$$p(x_k\mid y_{1:N},\theta ,f)$$

$$=\int p(x_k\mid x_{k+1},y_{1:N},\theta ,f)p(x_{k+1}\mid y_{1:N},\theta ,f)\mathrm{d}{x}_{k+1}$$

$$=\int p(x_k\mid x_{k+1},y_{1:k},\theta ,f)p(x_{k+1}\mid y_{1:N},\theta ,f)\mathrm{d}{x}_{k+1}.$$

$$p(x_k\mid x_{k+1},y_{1:k},\theta ,f)$$

$$=\frac{p(x_{k+1}\mid x_k,\theta ,f)p(x_k\mid y_{1:k},\theta ,f)}{p(x_{k+1}\mid y_{1:k},\theta ,f)}.$$

$$p(x_k\mid y_{1:N},\theta ,f)=p(x_k\mid y_{1:k},\theta ,f)\int \frac{p(x_{k+1}\mid x_k,\theta ,f)p(x_{k+1}\mid y_{1:N},\theta ,f)}{p(x_{k+1}\mid y_{1:k},\theta ,f)}\mathrm{d}{x}_{k+1}.$$

$$p(x_k\mid y_{1:k-1},\theta ,f)=\mathcal{N}\left(x_k\mid {\psi }^f,C_{\psi \psi }^f\right),$$

$$p(y_k\mid x_k)=\mathcal{N}\left(y_k\mid M{x}_k,C_{ϵϵ}\right),$$

$$p(x_k\mid y_{1:k},\theta ,f)=\mathcal{N}\left(x_k\mid {\psi }^a,C_{\psi \psi }^a\right),$$

$${\psi }^a={\psi }^f+{\left(M{C}_{\psi \psi }^f\right)}^T{\left[C_{ϵϵ}+M{C}_{\psi \psi }^f{M}^T\right]}^{-1}\left(y_k-M{\psi }^f\right),$$

$$C_{\psi \psi }^a=C_{\psi \psi }^f-{\left(M{C}_{\psi \psi }^f\right)}^T{\left[C_{ϵϵ}+M{C}_{\psi \psi }^f{M}^T\right]}^{-1}\left(M{C}_{\psi \psi }^f\right),$$

$$\overline{\psi }\approx \frac{1}{n}\sum _{j=1}^n{\psi }^j,\Psi =\left(\begin{array}{c}{\psi }^1-\overline{\psi },{\psi }^2-\overline{\psi },\cdots ,{\psi }^n-\overline{\psi }\end{array}\right),C_{\psi \psi }\approx \frac{1}{n-1}\Psi {\Psi }^T.$$

$${\left({\psi }^a\right)}^j=\overline{{\psi }^a}+{\left({\Psi }^a\right)}_j,$$

$$\overline{{\psi }^a}=\overline{{\psi }^f}+{\Psi }^f{\left(M{\Psi }^f\right)}^T{\left[(n-1)C_{ϵϵ}+M{\Psi }^f{\left(M{\Psi }^f\right)}^T\right]}^{-1}\left(y_k-M\overline{{\psi }^f}\right),$$

$${\Psi }^a={\Psi }^{f}V{\left[\mathbb{I}-\Sigma \right]}^{\frac{1}{2}}{V}^T,V\Sigma V^T={\left(M{\Psi }^f\right)}^T{\left[(N-1)C_{ϵϵ}+M{\Psi }^f{\left(M{\Psi }^f\right)}^T\right]}^{-1}M{\Psi }^f,$$

$$p(x_k\mid y_{1:k})=\mathcal{N}\left(x_k\mid {\psi }^s,C_{\psi \psi }^s\right),$$

$$\psi^{s}=\psi^{a}+\left(\big{(}C_{\psi\psi}^{f}\big{)}_{k+1}^{-1}GC_{\psi\psi}^{a}\right)^{T}\left[\psi_{k+1}^{s}-\psi_{k+1}^{f}\right]\quad,$$

$$C_{\psi\psi}^{s}=C_{\psi\psi}^{a}-\left(\big{(}C_{\psi\psi}^{f}\big{)}_{k+1}^{-1}GC_{\psi\psi}^{a}\right)^{T}\left[\big{(}C_{\psi\psi}^{f}\big{)}_{k+1}-\big{(}C_{\psi\psi}^{s}\big{)}_{k+1}\right]\left(\big{(}C_{\psi\psi}^{f}\big{)}_{k+1}^{-1}GC_{\psi\psi}^{a}\right)\quad.$$

$$\left(\psi^{s}\right)^{j}=\left(\psi^{a}\right)^{j}+\Psi^{a}\big{(}\Psi^{f}\big{)}_{k+1}^{-1}\left[\left(\psi^{s}\right)_{k+1}^{j}-\left(\psi^{f}\right)_{k+1}^{j}\right]\quad.$$

$$\frac{\mathrm{d}r}{\mathrm{d}t}=u(r),$$

Full text

††journal: Journal of Computational Physics

1 Introduction

A number of problems in computational physics involve the motion of interfaces, e.g. semiconductor manufacturing, multi-phase flows, crystal growth, groundwater flow, computer vision, grid generation and seismology [1]. In general, there are two approaches to calculating the motion of an interface [2]: In front-tracking methods, the interface is parameterized and discretized so that one follows the motion of the whole interface by tracking a sufficient number of points on the interface. In front-capturing methods, the interface is embedded into a function defined over the whole domain. One example of front-capturing methods are the so-called level-set methods, where the interface is embedded into a strictly monotonic function [3]. Both front-tracking and front-capturing methods have well-understood advantages, and hybrid methods exist to mitigate their respective disadvantages [4]. Nevertheless, it is worth mentioning that level-set methods provide a natural formulation for calculating the motion of an interface [1]: Level-set methods deal well with non-smooth features such as corners and cusps as well as topological merging and break-up. Furthermore, level-set methods are easily extended from two to three and higher dimensions. For more information on level-set methods, the reader is referred to many excellent expositions in the literature [5, 6, 1, 4, 2].

Despite the elegance of explaining physical phenomena by the motion of interfaces and calculating the motion of the interfaces by level-set methods, one has to remain aware that the assumption that manifolds are infinitely thin is often an asymptotic assumption. This is particularly true in fluid mechanics, which is governed by conservation laws [7]. One example is the kinematics of premixed flames [8]: Depending on the combustion regime, it may be assumed that a thin reactive-diffusive layer separates the burnt and unburnt gases. While the laminar flame speed, at which the premixed flame propagates from the burnt into the unburnt region, is a well-defined thermo-chemical property of the fuel-air mixture in a one-dimensional flow, it more generally depends on the balance between heat conduction and mass diffusion, whose effects vary with flame stretch and curvature. Moreover, the turbulent flame speed scales differently depending on the interactions between turbulence and combustion length scales. Despite the potential quantitative shortcomings in describing the kinematics of a premixed flame as the motion of an interface, it has been shown that this model successfully explains the linear and nonlinear dynamics observed in thermoacoustic instabilities of ducted premixed flames, which are relevant to the design of combustion chambers in jet and rocket engines [9, 10, 11, 12]. The discrepancies which arise when comparing to more faithful simulations or experiments are usually attributed to the unpredictable nature of turbulent flow or to uncertainties in the model and its parameters. Thus, it is relevant to assess the ability of a qualitative, physics-informed model to make quantitative, time-accurate predictions. The aim of this study is to develop a data-driven, statistically rigorous framework for state and parameter estimation in models using level-set methods, which has not been done before, and to apply it to an oscillating flame.

Inference over interfaces based on level-set methods has been performed, e.g. in the context of shape optimization, either by directly taking the functional derivative of the objective functional with respect to the level-set function [13, 14], or, equivalently, by embedding the level-set method into shape calculus [15, 16, 17]. As an alternative to these variational approaches, the framework for state and parameter estimation in this study is based on data assimilation [18]. Data assimilation finds the statistically optimal combination of model predictions and observations. It combines concepts from control theory, probability theory and dynamic programming [19, 20, 21]. The data assimilation technique used in this study is the ensemble Kalman filter [22, 23]. In the ensemble Kalman filter, a Monte-Carlo approach is used to represent the necessary statistics at every timestep [24], which makes it a computationally efficient technique in terms of storage requirements. Compared to other data assimilation techniques based on the Kalman filter, e.g. the extended Kalman filter [25, 26], the ensemble Kalman filter is found to be particularly robust with respect to larger nonlinearities [18]. This is relevant for level-set methods, due to strongly nonlinear events such as cusp formation, topological merging and break-up. A practical advantage of the ensemble Kalman filter is its non-intrusive implementation with little effort required for its parallelization. The ensemble Kalman filter has been successfully applied to a number of problems in fluid mechanics: turbulent near-wall flow [27]; transonic flows around airfoils and wings [28]; viscous flow around a cylinder [29]; model uncertainties in Reynolds-averaged Navier-Stokes (RANS) equations [30]; vortex models of separated flow [31]; and extinction and reignition dynamics in turbulent non-premixed combustion [32].

In this study, the ensemble Kalman filter is combined with a narrow-band level-set method and a fast marching method to form a computationally efficient level-set data assimilation framework [40, 41]. In analogy to the distinction between front-tracking and front-capturing methods in describing the motion of an interface, data can also be assimilated according to various paradigms, e.g. as demonstrated by [42], [43] and [44]. For the various paradigms, it is not a priori clear which will yield superior results. For this reason, the theory and the derivations behind our level-set data assimilation framework are worked out in detail. The paper is structured as follows: In Section 2, the ensemble Kalman filter and the ensemble Kalman smoother are derived within the context of Bayesian inference. The estimation problems in data assimilation and parameter estimation are then formulated in terms of probability distributions. In Section 3, various formulations for the laws of motion are discussed, including the Hamilton-Jacobi equation. It is shown that the solutions to the Hamilton-Jacobi equation, the so-called generating functions, form a natural state space for data assimilation. The level-set data assimilation framework based on the Hamilton-Jacobi equation is verified in one-dimensional and two-dimensional examples. In Section 4, statistical inference using the level-set data assimilation framework is demonstrated on the nonlinear dynamics of a ducted premixed flame. The results and insights of this study are summarized in the conclusions (Section 5).

2 Data assimilation and parameter estimation

Data assimilation and parameter estimation are here treated as problems in statistical inference. Statistical inference quantifies the degree of belief (or confidence) in a physical model, the parameters that it receives and the states that it predicts. Statistical inference follows a probabilistic formulation, which provides precise definitions for the different tasks addressed in inference, which include filtering, smoothing and prediction. Probabilistic formulations for statistical inference include frequentist inference and Bayesian inference. In frequentist inference, an error functional is defined to measure the statistical distance between a candidate solution and the available data. Inference becomes an optimization problem to minimize this error functional. In Bayesian inference, existing knowledge is quantified in the form of a probability distribution over candidate solutions. When data becomes available, the probability distribution is updated, effectively combining the existing knowledge with the data. For normal distributions under linear dynamics, both formulations of statistical inference give equivalent results [18]. Under these circumstances, frequentist inference may be considered more accessible because it allows the application of familiar tools from convex optimization [45]. Nevertheless, Cox’s axioms demonstrate that probability theory, as used in Bayesian inference, gives a natural formulation for inference in the general case [46]. All definitions and derivations in this section are given in terms of Bayesian inference (A). For more details, the reader may refer to [18] or [47].

2.1 Probabilistic state space model

The state of a system at a timestep $k$ is uniquely defined by the state vector $x_{k}$. The evolution of the state is governed by a physical model $f$ and its parameters $\theta$. Thus, our belief in the state $x_{k}$ depends exclusively on our belief in (i) the previous state $x_{k-1}$ as well as our choice of (ii) the physical model $f$ and (iii) its parameters $\theta$. At the same time, the state vector $x_{k}$ is compared to noisy observations $y_{k}$ through a measurement operator $M$. Formulated in terms of a probabilistic state space model, we obtain

[TABLE]

[TABLE]

In brief, the transition from one state to the next is governed by the operator $G$. The operator $G$ is derived from the physical model $f$, and depends on the parameters $\theta$. The states $x_{k}$ and the observations $y_{k}$ are considered realizations inside their respective probabilistic state spaces. The degrees of belief in each are subject to the conditional probability distributions $p(x_{k}\mid x_{k-1},\theta,f)$ and $p(y_{k}\mid x_{k})$ respectively. The degree of belief in the state $x_{k}$ is conditional on the previous state $x_{k-1}$ and the choices in the physical model $f$ and its parameters $\theta$. The degree of belief in the observation $y_{k}$ is conditional on the true state $x_{k}$.

The rules of the probabilistic state space model may be formalized as follows [47]. Firstly, we assume that the physical model may be described by a Markov chain [24]. This means that the belief in a state depends only on the belief in the previous state:

The future is independent of the past given the present:

[TABLE]

where $x_{0:k-1}$ denotes the states at all timesteps from 0 to $k-1$, and $y_{1:k-1}$ denotes the observations at all timesteps from 1 to $k-1$. 2. 2.

The past is independent of the future given the present:

[TABLE]

where $x_{k+1:N}$ and $y_{k+1:N}$ respectively denote the states and observations at all timesteps from $k+1$ to the final timestep $N$.

Secondly, observations are assumed to be conditionally independent in time. The probability of an observation depends only on the current state:

[TABLE]

In Fig. 1, the relationship between states and observations is shown as well as the roles of models, parameters and measurement.

The goal of data assimilation is to find the joint probability distribution $p(x_{0:N},y_{1:N},\theta,f)$. The probabilistic state space spans all states $x_{k}$ from timestep [math] to $N$ and all observations $y_{k}$ from timestep $1$ to $N$, as well as the physical model $f$ and its parameters $\theta$. This joint probability distribution gives a complete statistical description. In principle, all probability distributions of interest may be derived from this joint probability distribution (Table 1). In practice, it is difficult to compute this probability distribution because of the high dimensionality of its probabilistic state space spanning multiple timesteps [48]. Therefore, state estimation focuses on the more direct computation of conditional probability distributions over a single timestep (Table 2).

In the probabilistic formulation, data assimilation is easily extended to account for parameters. In combined state and parameter estimation, the state is augmented by the parameters so that they become subject to the same inference (Eq. (1)):

[TABLE]

The tasks in combined state and parameter estimation are given in Table 3. The filtered and smoothed distributions in the parameters $\theta_{k}$, $p(\theta_{k}\mid y_{1:k},f)$ and $p(\theta_{k}\mid y_{1:N},f)$ respectively, are retrieved by marginalizing the states $x_{k}$ (A). Note that the parameters $\theta_{k}$ are now time-dependent as the system traverses different regimes in state space. This turns the strongly constrained parameter estimation into a weakly constrained combined state and parameter estimation [18]. Thus, the results of marginalizing the probability distributions in combined state and parameter estimation (Table 3) are not strictly equivalent to the solutions of parameter estimation (Table 1).

2.2 Bayesian filtering and smoothing

For the filtering problem, Bayes’ rule gives

[TABLE]

The prediction $p(x_{k}\mid y_{1:k-1},\theta,f)$ is given by the Chapman-Kolmogorov equation:

[TABLE]

The Chapman-Kolmogorov equation requires the inverse from the previous timestep $k-1$. In general, it is solved either numerically [49], analytically (Theorem 3) or by a Monte-Carlo simulation (Theorem 4). The steps in the Bayesian filter may be summarized as follows:

Theorem 1** (Bayesian filter).**

**

Prediction step:

[TABLE] 2. 2.

Update step:

[TABLE]

At the timestep $k=N$, the filtered and smoothed distributions are identical (Tab. 2). If the smoothed distribution at a timestep $k+1$ is known, the smoothed distribution at the previous timestep $k$ is also known due to the Markov chain properties of the probabilistic state space model (Eq. (3), (4)). This may be formalized as follows:

[TABLE]

$p(x_{k+1}\mid y_{1:N},\theta,f)$ is the smoothed distribution from the subsequent timestep $k+1$. $p(x_{k}\mid x_{k+1},y_{1:k},\theta,f)$ is computed via Bayes’ rule:

[TABLE]

Note that it involves the filtered distribution from the current timestep $k$. The steps in the Bayesian smoother may be summarized as follows:

Theorem 2** (Bayesian smoother).**

**

Forward sweep: Bayesian filter (Theorem 1). 2. 2.

Backward step:

[TABLE]

Note that both the Bayesian filter and smoother are sequential in nature. In the Bayesian filter, the predicted distribution $p(x_{k}\mid y_{1:k-1},\theta,f)$ at a timestep $k$ (Eq. (12)) mainly depends on the filtered distribution $p(x_{k-1}\mid y_{1:k-1},\theta,f)$ from the previous timestep $k-1$ (Eq. (13)). In the Bayesian smoother, the smoothed distribution $p(x_{k}\mid y_{1:N},\theta,f)$ at a timestep $k$ mainly depends on the smoothed distribution $p(x_{k+1}\mid y_{1:N},\theta,f)$ from the subsequent timestep $k+1$ (Eq. (19)). The existence of sequential algorithms for the computation of filtered and smoothed distributions significantly reduces the complexity of data assimilation [48].

2.3 The Kalman filter and the Rauch-Tung-Striebel smoother

Two additional assumptions are introduced to make the computation of filtered and smoothed distributions feasible. Firstly, the prior and the likelihood in the update step of the Bayesian filter are assumed to be normal (Eq. (13)):

[TABLE]

[TABLE]

where $\mathcal{N}$ denotes a normal distribution with respective mean and covariance matrix. The mean of the prior is denoted by $\psi^{f}$, its covariance matrix by $C_{\psi\psi}^{f}$, and the covariance matrix of the likelihood, also known as the observation error, by $C_{\epsilon\epsilon}$. From Eq. (20) and (21), it follows that the filtered distribution $p(x_{k}\mid y_{1:k},\theta,f)$ is normal (Eq. (13)). Secondly, the operator $G$ (Eq. (1)) in the prediction step (Eq. (12)) is assumed to be linear in $x$. The result is the well-known Kalman filter [50, 51]:

Theorem 3** (Kalman filter).**

[TABLE]

[TABLE]

[TABLE]

where the superscript $f$ denotes ’forecast’ (everything pertaining to the prediction), and the superscript $a$ denotes ’analysis’ (everything pertaining to the update). If the physical model is nonlinear, the prediction of the covariance matrix $C_{\psi\psi}^{f}$ for the Kalman filter may be generalized in several ways. In the extended Kalman filter, the predicted covariance matrix $C_{\psi\psi}^{f}$ is computed by linearizing $f$ [20]. In strongly nonlinear dynamical systems, the predictions are found to be poor [18]. Higher-order extended Kalman filters are available [25]. Nevertheless, drawbacks include their significant storage requirements, which increase exponentially with the order of approximation. An alternative is the ensemble Kalman filter [22, 23]. Instead of a mean $\psi$ and a covariance matrix $C_{\psi\psi}$, a distribution is represented by a sample $\left(\psi^{j}\right)_{j=1\dots n}$. During the prediction step, the ensemble members $\psi^{j}$ evolve in time independently. Before the update step, the statistics may be recovered from the sample as follows:

[TABLE]

The sample covariance matrix $C_{\psi\psi}$ involves division by $n-1$ instead of $n$ in order to avoid a sample bias.

Various implementations of the ensemble Kalman filter exist, which differ in the update step. In the straightforward implementation of the ensemble Kalman filter [22], every ensemble member is individually updated (Eq. (23)). It can be shown that the observations must be randomly perturbed in order to guarantee a statistically consistent analysis scheme [23]. In order to avoid the introduction of randomly generated numbers, the square-root filter is used here [52]. The square-root filter belongs to a larger family of ensemble Kalman filters called ensemble-transform Kalman filters [53, 54]. Unlike the straightforward implementation of the ensemble Kalman filter, the mean and the deviations of the ensemble members are updated. This requires the singular value decomposition of a symmetric, positive, semi-definite matrix ($V\Sigma V^{T}$, where $V$ is orthonormal and $\Sigma$ diagonal), but no spurious errors due to the random perturbation of the observations are introduced.

Theorem 4** (Square-root filter).**

[TABLE]

[TABLE]

[TABLE]

where $\mathbb{I}$ is the identity matrix.

Following Eq. (20) and (21), the Bayesian smoother (Theorem 2) becomes the Rauch-Tung-Striebel smoother, also known as Kalman smoother (Theorem 5). It can be shown that the smoothed distribution becomes a normal distribution with mean $\psi^{s}$ and covariance matrix $C_{\psi\psi}^{s}$. In order to again avoid the shortcomings of assuming linearity, an ensemble Kalman smoother is presented here (Theorem 6).

Theorem 5** (Rauch-Tung-Striebel smoother).**

[TABLE]

[TABLE]

[TABLE]

Theorem 6** (Ensemble Kalman smoother).**

[TABLE]

3 Level-set methods

When data assimilation is applied to the motion of an interface, the question arises as to what constitutes its probabilistic state space. Refining the distinction between front-tracking and front-capturing methods, there are at least three ways to view the motion of an interface [1, 4]:

Geometric view.

The interface is parameterized and discretized so that one follows the motion of the whole interface by solving the laws of motion for a sufficient number of points on the interface.

Set-theoretic view.

A characteristic function is defined over the whole domain. The characteristic function assumes one of two values, depending on whether the point at the location in question is inside or outside the region enclosed by the interface.

Analytic view (also known as “analysis view” [4]).

A level-set function is defined over the whole domain. The interface is reconstructed by identifying the position of a particular level set.

Beyond the description of motion, these three views also pertain to the assimilation of data. In fact, all three views have been employed in earlier level-set data assimilation frameworks: While using a level-set method when solving for the motion of the interface, [43] take each entry in the innovation vector $y_{k}-M\psi^{f}$ (Eq. (22)) as the distance between one point on the predicted interface and one point on the observed interface. This framework corresponds to the geometric view. [44] compute $y_{k}$ and $\psi^{f}$ in the innovation vector from a progress variable. The progress variable assumes values between zero and one, with most intermediate values assumed near the interface. This framework is an approximation of the set-theoretic view. [42] compute $y_{k}$ and $\psi^{f}$ from a level-set function. This framework corresponds to the analytic view, and comes closest to the level-set data assimilation framework presented in this study. The motion of an interface depends on the values of the level-set function in the immediate environment of the interface, whereas the values of the level-set function are not unique away from the interface. The level-set data assimilation framework presented in this study (i) introduces an additional constraint in agreement with the the analytic view, (ii) validates the choice of the constraint for two analytical test cases in one and two dimensions respectively, and (iii) addresses shortcomings in the aforementioned frameworks.

In order to construct the probabilistic state space, the derivation of the level-set method is revisited in this section. The centerpiece of this derivation is the Hamilton-Jacobi equation. In theory, the level-set method may alternatively be derived as the transport of a passive scalar quantity [8, Chapter 2]. Although the value of the level-set function is well defined at the interface, the choice of values for the level-set function away from the interface is in general not unique, and would require an ad-hoc constraint. It is demonstrated how the necessary constraint to the Hamilton-Jacobi equation naturally follows from the choice of phase space. Finally, it is shown that the solutions to the Hamilton-Jacobi equation, the so-called generating functions, form the appropriate probabilistic state space for data assimilation. Two complimentary level-set algorithms are combined to solve the Hamilton-Jacobi equation computationally efficiently: the narrow-band method [40] and the fast-marching method [41]. Along with the ensemble Kalman filter (Theorem 4), the three algorithms form the backbone of our level-set data assimilation framework.

3.1 Hamilton-Jacobi equation

The laws of motion of an interface shall be given by

[TABLE]

where $r$ is the position of one point on the interface, and $u$ is the velocity. In Hamiltonian mechanics, motion is described in phase space using generalized coordinates and momenta [55, Chapter 3]. With interfaces in mind, a natural choice for the phase space is to use the position $r$ as generalized coordinates and the normal vector $n$ as generalized momenta.

For the Hamilton-Jacobi equation, the generating function $G(r,t)$ shall satisfy the following properties:

The initial interface $r(0)$ is given by a level set $G(r(0),0)=\mathrm{const}$.

• 2.

The moving interface is identified as the level set $G(r(t),t)=G(r(0),0)$.

• 3.

Huygens’ principle [55, Chapter 9] states that $\nabla G=n$. Given the choice of phase space, this is a necessary relationship between the generating function and the generalized coordinates and momenta.

The first two properties are present in every front-capturing method [8]. The third property, Huygens’ principle, follows from Hamiltonian mechanics (B). Physically speaking, every point on the interface is interpreted as the source of a wave. As the waves propagate through the domain, they connect each point on the interface to every other point in the domain. Geometrically speaking, Huygens’ principle gives an eikonal field [56].

It remains to show how the generating function translates into a level-set method. The Lagrangian $\mathcal{L}(r,\dot{r},t)$ is given by (B, Lemma 1)

[TABLE]

[TABLE]

[TABLE]

The Hamiltonian $\mathcal{H}(r,n,t)$, subject to the same constraints as the Lagrangian $\mathcal{L}(r,\dot{r},t)$ (Eq. (35), (36)), is given by (B, Lemma 2)

[TABLE]

Finally, the Hamilton-Jacobi equation is given by (B, Theorem 9)

[TABLE]

The solution to the Hamilton-Jacobi equation, the generating function $G(r,t)$, is defined over the whole domain. As such, it represents the state of the interface in the probabilistic state space. Eq. (38) is solved near the interface. Away from the interface, the generating function is constrained by (Eq. (36))

[TABLE]

Eq. (39) is an eikonal equation [56]. The solutions to this eikonal equation are signed distance functions.

3.2 Level-set data assimilation framework

The level-set data assimilation framework presented in this study combines three algorithms: (i) the narrow-band method to solve the Hamilton-Jacobi equation near the interface [40], (ii) the fast-marching method to extend the eikonal field from the narrow band to the whole domain [41], and (iii) the ensemble Kalman filter and smoother to assimilate data [18].

The narrow-band method forms the backbone of the level-set data assimilation framework [40]. It solves the Hamilton-Jacobi equation (Eq. 38) in a narrow band near the interface. Using a narrow band reduces the computational cost by one order of magnitude compared to finding the generating function over the whole domain [57]. In a numerical simulation, the spatial resolution shall be denoted by $\Delta x$ and the temporal resolution by $\Delta t$. The temporal resolution has to satisfy the Courant-Friedrich-Lewy condition $u\Delta t<\Delta x$ at every point on the interface [40]. The narrow-band method consists of the following steps (Fig. 2, left):

Initialize.

Before the first timestep, the signed distance field has to be known in a sufficiently wide neighborhood of the initial interface. It may be obtained from either a previously computed generating function or exact knowledge of the position of the initial interface.

Localize.

At the timestep $k$, a band $T^{k}=\{r:G(r,t_{k})<\gamma\Delta x\}$ is formed, where $\gamma$ is a constant chosen according to the width of the narrow band. A second, enveloping band $N^{k}=\{r:G(r+r^{\prime},t_{k})<\gamma\Delta x\,,\,\exists|r^{\prime}|<\Delta x\}$ is formed. It follows from the Courant-Friedrich-Lewy condition that the band $N^{k}$ contains the band $T^{k+1}$ formed at the next timestep $k+1$.

Advance.

The generating function over the band $N^{k}$ is evolved by one timestep. Note that the solution $\widetilde{G}^{k+1}(r)$ is not the generating function at the next timestep $k+1$. While the position of the interface correctly coincides with the level set of the solution $\widetilde{G}^{k+1}(r)$, Huygens’ principle is no longer satisfied [58].

Reinitialize.

The following Hamilton-Jacobi equation is solved over the band $N^{k}$ until steady state is reached [58]:

[TABLE]

[TABLE]

The steady-state solution $\widetilde{G}(r,\tilde{t}\to\infty)$ gives the generating function $G(r,t_{k+1})$ over the band $N^{k}$ at the timestep $k+1$. Note that this Hamilton-Jacobi equation is a numerical continuation of the eikonal equation (Eq. (39)), where $\widetilde{t}$ denotes a pseudo-time [59].

When observations are available, the prediction step is followed by an update step (Fig. 2, right). Like the prediction step, the update step involves localization and reinitialization so that the narrow band follows the position of the interface. Instead of the advancement step, there is an extension step and an assimilation step:

Extend.

While a narrow band is sufficient to capture the motion of an interface, data assimilation has to be performed in one probabilistic state space for all interfaces. In a numerical simulation, the state comprises the values of the generating function at every grid point. The fast-marching method extends the generating function from the narrow band to the whole domain, thus returning a well-defined state. The complexity of the fast-marching method is $\mathcal{O}(m^{2}\log(m))$ in two dimensions or $\mathcal{O}(m^{3}\log(m))$ in three dimensions, where $m$ denotes the number of grid points in one dimension [41]. While the generating function over the whole domain is relatively expensive to compute, the fast-marching method is only called when observations are available. If combined with the narrow-band method (complexity $\mathcal{O}(m)$ or $\mathcal{O}(m^{2})$ [40]), this level-set method becomes affordable overall.

Assimilate.

The ensemble Kalman filter (Theorem 4) and/or smoother (Theorem 6) are applied. The assimilation step is followed by a localization step and a reinitialization step in order to prepare the narrow band for the next prediction step. Note that the result of data assimilation is in general not an eikonal field. While the band localized after data assimilation may in theory be too narrow, it can be shown that the slopes in the filtered solutions are in expectation less than or equal to unity (C).

Combining the building blocks of the level-set method and data assimilation gives the level-set filtering and smoothing frameworks (Fig. 3). In the following subsections, the update step in the level-set filtering framework is verified using two analytical test cases: In a one-dimensional test case, the level-set filtering framework, which represents the analytic view of data assimilation, is contrasted with the geometric and set-theoretic views of data assimilation. It is shown that the analytic view of data assimilation works correctly, which sets the stage for its generalization to higher dimensions. Afterwards, the level-set filtering framework is applied to a two-dimensional test case, and the effect of the observations on the shape of the interface is discussed. It is shown how various parameters such as the number of observations and the observation error may affect the performance of the level-set data assimilation framework.

3.3 One-dimensional test case

The predicted and observed positions $x$ and $y$ of an interface in one dimension are respectively given by (Fig. 4)

[TABLE]

[TABLE]

In one dimension, an interface reduces to a point. The normal vector on the interface is unique up to its orientation. Geometric data assimilation using the Kalman filter (Theorem 3) reduces to optimal interpolation of the positions $\mu_{\psi}$ and $\mu_{\epsilon}$ (Eq. (23)):

[TABLE]

The filtered position $x$ is closer to the observation than to the prediction in accordance with the variances $C_{\psi\psi}$ and $C_{\epsilon\epsilon}$ (Fig. 4, left).

For set-theoretic data assimilation, the characteristic function $\chi$ is assumed to be bounded by $\pm 1$ with the interface at $\chi=0$:

[TABLE]

The characteristic function $\chi$ follows a Bernoulli distribution, which is the appropriate probability distribution for a random variable with two possible outcomes, at every location $x$ [60, Chapter 6]:

[TABLE]

where $\Phi$ denotes the cumulative density function of the normal distribution. The characteristic function $\chi$ does not follow a normal distribution whereas the Kalman filter requires the mean $\overline{\chi}$ and the covariance function $\mathcal{C}_{\chi\chi}$:

[TABLE]

Note that the mean $\overline{\chi}$ is not a characteristic function. Set-theoretic data assimilation using the Kalman filter (Theorem 3) gives

[TABLE]

The interface is localized at $x\approx 0.139$ (Fig. 4, middle). The filtered position $x$ is closer to the prediction than to the observation although the prediction variance $C_{\psi\psi}$ is larger than the observation variance $C_{\epsilon\epsilon}$. Note that the filtered position $x$ is independent of the variance $C_{\epsilon\epsilon}$.

For analytic data assimilation, the predicted and observed positions $x$ and $y$ are viewed as points on the level sets of the predicted and observed generating functions respectively. It follows from Eq. (42) that the predicted generating function is given by

[TABLE]

The mean $\overline{G}$ and the covariance function $\mathcal{C}_{GG}$ are given by

[TABLE]

Analytic data assimilation using the Kalman filter (Theorem 3) on the generating function gives

[TABLE]

The filtered position $x$ is given by the level set of the filtered generating function $G$ (Eq. (52)), which is identical to the result from geometric data assimilation (Fig. 4, right).

In summary, geometric, set-theoretic and analytic data assimilation are compared. Set-theoretic data assimilation gives implausible results in one dimension, and is thus discarded from any further consideration. Geometric and analytic data assimilation give consistent results in one dimension. Nevertheless, the generalization of geometric data assimilation to higher dimensions is not straightforward because it requires a one-to-one correspondence between the predicted interface and the observations. In the next test case, the generalization of analytic data assimilation, as the most feasible framework, to higher dimensions is investigated.

3.4 Two-dimensional test case

The distance field $S$ of a corner in $(x_{0},y_{0})$ is given by

[TABLE]

For prediction, the coordinates of the corner $(x_{0},y_{0})$ are taken to be independent and normally distributed:

[TABLE]

A sketch of the corner and its median distance field $S(x,y;0,0)$ are given in Fig. 5.

For analytic data assimilation, the predicted position $x$ is again viewed as the level set of a predicted generating function. It follows from Eq. (53) and (54) that the predicted generating function is given by

[TABLE]

where the mean $\overline{G}$ and the covariance function $C_{GG}$ are given by

[TABLE]

where $\mathbb{E}_{x_{0},y_{0}}$ denotes the expected value over the random variables $x_{0}$ and $y_{0}$, and the ’$\times$’ symbol denotes the Cartesian product of two random variables. The mean $\overline{G}$ and the covariance function $C_{GG}$ are computed from Eq. (57) and (60) using Gauss-Hermite quadrature [61]. Alternatively, the mean and the variance function can be computed from the following one-dimensional integrals on the unit circle:

[TABLE]

The auxiliary functions $g$ and $h$ are given by

[TABLE]

where $\varphi$ and $\Phi$ denote the probability density function and cumulative distribution function of the normal distribution. The natural coordinates $(v,w)$ are given by

[TABLE]

The mean and the standard deviation of the predicted generating function are shown in Fig. 6. Artefacts from the interaction between the nonlinearity of distance fields and the assumption of normal distributions are visible: The mean of the predicted generating function (Fig. 6, left) is rounded in the corner whereas distance fields are supposed to have right angles (Fig. 5). The standard deviation of the predicted generating function (Fig. 6, right) assumes non-unity values along $\theta=45^{\circ}$ although the coordinates of the corner are independent and normally distributed (Eq. (54)). The effects of the artefacts in the predicted generating function are now examined in various observation scenarios.

Unlike the one-dimensional example, where one point fully characterizes an interface, observations in two (and higher) dimensions vary in resolution, ranging from the position of one point on the interface to a full sampling of the interface. In Fig. 7, the influence of one observation depending on its location is investigated. It is evident that one observation is insufficient to represent a corner. Nevertheless, analytic data assimilation gives qualitatively meaningful results in that the observation shifts the vertical and horizontal edges proportionally depending on the location of the observation. This is fundamentally different from geometric data assimilation (Fig. 4, left), where every point on the interface is required to correspond to an observation. This makes the presented (analytic) level-set data assimilation framework fully front-capturing (rather than front-tracking), both with respect to the description of the motion of the interface and the assimilation of data. As the number of observation points increases, the right angle in the corner is retrieved (Fig. 8).

It is worth pointing out the probabilistic nature of the presented level-set data assimilation framework. Similar to soft clustering [60], the points on the predicted interface are not associated with distinct observations but with all observations to varying degrees. If we consider $\psi^{a}$ to be a function of $\psi^{f}$, $C_{\psi\psi}^{f}$, $y$, $C_{\epsilon\epsilon}$ and $M$ (Eq. (23)), the sensitivity to data is given by

[TABLE]

where the columns represent the sensitivities to individual observations. The shift in the zero-level set (Fig. 8) is the superposition of the contributions due to the individual observations. For example, while the points on the vertical edge far away from the vertex are mainly affected by the uppermost observation (Fig. 9(a)), the other observations have a vanishing effect on them as well (Fig. 99(b)-9(d)). Nevertheless, the shape deformations are more localized compared to Fig. 7. While Eq. (65) does not formally depend on the observations $y$, it includes information about the number and the locations of the observations signified by the measurement operator $M$. Note that each observation of interest leaves the values of the generating function at the other observation locations intact.

In summary, analytic data assimilation is applied to a prototypical two-dimensional shape, a corner. Within the level-set data assimilation framework, the differences between distance fields, predicted generating functions and filtered generating functions are illustrated. In particular, the effects of the number of observation points and their locations are studied. In the following section, the level-set data assimilation framework is applied to a more realistic experiment. In general, the predicted generating function is taken as a composition of straight sections, which exhibit the behavior discussed for the one-dimensional test case, and strongly curved sections insufficiently resolved by the observations, which exhibit the behavior discussed for the two-dimensional test case. The theoretical insights from this section will be referenced throughout the next section.

4 Twin experiments for a premixed flame inside a duct

Thermoacoustic instabilities are a persistent challenge in the design of jet and rocket engines: Velocity and pressure oscillations inside the combustion chamber interact with the flame and cause an unsteady heat release rate. If moments of higher heat release rate coincide with moments of higher pressure (and lower heat release rate with lower pressure), acoustic oscillations arise. This can lead to large-amplitude oscillations causing structural damage in the jet or rocket engine [62, 63].

The so-called $G$-equation model is a reduced-order model to study the flame dynamics which lead to heat release rate perturbations [9, 10]. The premixed flame is modeled as an interface captured by a level set. The velocity of a point in the level set is the sum of the flame speed, which is normal to the interface and points towards the unburnt gas, and the underlying flow field. The underlying flow field includes both hydrodynamic and acoustic contributions. The $G$-equation model is used to compute flame transfer functions (FTF) or flame describing functions (FDF), which give the heat release perturbation to a given velocity perturbation. Coupled with linear acoustics models, the $G$-equation model has been very successful in qualitatively characterizing the linear and nonlinear dynamics of self-excited thermoacoustic oscillations [11, 12, 64, 65].

We demonstrate our level-set data assimilation framework by performing twin experiments for the $G$-equation model applied to a ducted premixed Bunsen flame under acoustic forcing. In our twin experiments, both the model predictions and the observations come from $G$-equation simulations, but with different sets of parameters [66, Chapter 9]. This leads to uncertainties in the parameters and, as a result, in the states. In the absence of uncertainties in the model, the twin experiment is an important benchmark in the quantitative assessment of a data assimilation framework because it becomes possible to compare filtered and smoothed solutions to a reference solution at all times [66, Chapter 8]. An example for the application of our level-set data assimilation framework to more realistic observations from a direct numerical simulation can be found in the proceedings of the 2018 CTR Summer Program [67].

4.1 Reference solution

The $G$-equation is given by [8, Chapter 2]

[TABLE]

where the underlying flow field is denoted by $u(r,t)$ and the laminar flame speed by $s_{\mathrm{L}}$. For $|\nabla G|=1$, the $G$-equation is formally equivalent to the Hamilton-Jacobi equation (Eq. 38). The underlying flow field $u(r,t)$ is a superposition of the steady, inviscid flow field $U$ and the time-dependent perturbation velocity field $u^{\prime}(r,t)$. In the reduced-order model of $u^{\prime}(r,t)$, the perturbation has an amplitude $\varepsilon U$, and travels downstream at the phase speed $U/K$, where $K$ and $\varepsilon$ are model parameters [69]. In the axisymmetric case, the cylindrical components of the perturbation velocity field are given by

[TABLE]

where the wavenumber $k$ and the angular frequency $\omega$ satisfy the dispersion relation $\omega(k)=Uk/K$. The perturbation velocity field $u^{\prime}(r,t)$ is, mathematically speaking, divergence-free, and, physically speaking, satisfies the continuity equation.

In Fig. 10, a sketch of the ducted premixed flame is shown. The ducted premixed flame is forced at an angular frequency of $277.78\,\mathrm{rad}/\mathrm{s}$. Tab. 4 shows the measurements of the burner and the tube. The $G$-equation is numerically solved using the narrow-band level-set method with distance reinitialization [40, 58]: The computational domain is discretized using a weighted essentially non-oscillatory (WENO) scheme in space and a total-variation diminishing (TVD) version of the Runge-Kutta scheme in time, which give up to fifth-order accuracy in space and third-order accuracy in time [70, 71, 72]. At the base of the flame, a rotating boundary condition is used [73]. The $G$-equation solver has been verified in a number of studies [74, 11, 12, 75, 65]. These studies showed that the $G$-equation reproduces the dynamics of premixed flames only qualitatively, not quantitatively. With our level-set data assimilation framework, significant improvements on quantitative predictions are achieved. In axial and radial coordinates, a uniform $401\times 401$ Cartesian grid is used, which corresponds to $30\,\mu\mathrm{m}\times 180\,\mu\mathrm{m}$ grid cells. An advection Courant-Friedrich-Lewy (CFL) number of 0.02 is chosen, which corresponds to a timestep of 3.6 $\mu$s or approximately 6,283 timesteps over one period.

Snapshots of the reference solution to the $G$-equation for the ducted premixed flame are shown in Fig. 11. The flame is attached to the burner lip, while the perturbations travel from the base of the flame to the tip. When the perturbations are sufficiently large, a fuel-air pocket pinches off. In our reduced-order model, the perturbations are mainly governed by two non-dimensional parameters $K$ and $\varepsilon$, which govern phase speed and amplitude respectively [69]. In practice, neither parameter is accurately known a priori, which is a major source of uncertainty.

4.2 Combined state and parameter estimation

Three variations of the twin experiment are performed:

Monte-Carlo experiment.

An ensemble of 32 $G$-equation simulations is performed. Each simulation has a different set of parameters $K$ and $\varepsilon$. They are independently sampled from a normal distribution with a standard deviation of 20 %.

State estimation.

In addition to the procedure described for the Monte-Carlo experiment, the ensemble Kalman filter (Theorem 4) is applied every 1,000 timesteps to update the discretized generating function $G$ while leaving each set of parameters $K$ and $\varepsilon$ unaltered. The observations are extracted from the reference solution (Section 4.1): The measurement operator $M$ is an indicator matrix which identifies all grid points in the reference solution adjacent to its zero-level set. All aforementioned grid points are treated as observations of the flame surface. As such, all entries in the vector $y_{k}$ are set to 0. The observation error $\sigma_{\epsilon}$ is set to 60 $\mu$m, which corresponds to two grid cells.

Combined state and parameter estimation.

In addition to the procedure described for the state estimation, the discretized generating function $G$ is augmented by appending the parameters $K$ and $\varepsilon$ to the state vector (Eq. 6). Thus, both the state and the parameters are updated whenever data in the form of observations from the reference solution is assimilated.

In each variation of the twin experiment, the $k$-th entry in the state vector $\psi$ is marginally distributed according to

[TABLE]

The mean $\overline{\psi}[k]$ and the variance $C_{\psi\psi}[k,k]$ are computed from Eq. (25). Explicitly, the likelihood for the flame surface to be found at the location of the $k$-th entry is given by

[TABLE]

Alternatively, the logarithm of the normalized likelihood is given by [76]

[TABLE]

In Fig. 12, the logarithm of the normalized likelihood is shown for the three variations of the twin experiment. Its zero-level set gives the maximum-likelihood location of the flame surface. The more negative the value at a location is, the less likely the flame surface is to be found there. In Fig. 12(a), the logarithm of the normalized likelihood is shown for the Monto-Carlo experiment. As the perturbation travels from the base of the flame to the tip, the high-likelihood region for the location of the flame surface spreads out. The high-likelihood region is largest when fuel-air pockets pinch off, which represents maximal uncertainty. In Fig. 12(b), the logarithm of the normalized likelihood is shown for the twin experiment with state estimation. A qualitative comparison to Fig. 12(a) shows that the high-likelihood region for the location of the flame surface is significantly tighter. While the high-likelihood region still grows as the perturbation travels, the regular assimilation of data suppresses it to the vicinity of the observed flame surface. In Fig. 12(c), the logarithm of the normalized likelihood is shown for the twin experiment with combined state and parameter estimation. After few data assimilation cycles, the state-augmented ensemble Kalman filter has learnt the parameters $K$ and $\varepsilon$ almost exactly. Knowledge of the parameters enables highly precise predictions of the location of the flame surface, even during pinch-off events. Unlike in Fig. 12(a) and 12(b), no growth in the high-likelihood region is qualitatively discernable.

A global, more quantitative measure for the uncertainty in the location of the flame surface is the root mean square (RMS) error, which is defined as the square-root of the trace of the covariance matrix of the ensemble (Eq. (25)):

[TABLE]

In Fig. 13, the RMS error is plotted over time for the three variations of the twin experiment from $t=0$ for eight cycles. The RMS error is initially zero in all three twin experiments because all members in each ensemble share the same initial condition. In the Monte-Carlo experiment (Fig. 13, blue line), the RMS error subsequently grows until it reaches a high-uncertainty plateau. Momentary spikes in the uncertainty occur approximately every 6000 timesteps, and coincide with the pinch-off events observed in the reference solution (Fig. 11). With state estimation (Fig. 13, orange line), the assimilation of data regularly suppresses the uncertainty as qualitatively observed in Fig. 12(b). As a result, the predictions of the location of the flame surface significantly improve. With combined state and parameter estimation (Fig. 13, green line), the predictions further improve. Beginning with the first instance of data assimilation at timestep 1,000, the uncertainty steadily decreases until it reaches a low-uncertainty plateau after approximately 10,000 timesteps. At this point, the state-augmented ensemble Kalman filter has optimally calibrated the parameters, and the state has reached a statistically stationary state.

To gain insight into the effect of combined state and parameter estimation, it is instructive to analyze the probability distributions in the parameters. In Fig. 14, the marginal distributions in the normalized residuals of $K$ and $\varepsilon$ are shown, signified by their means and three-sigma confidence levels at every timestep. By timestep 15,000, no improvement is visible in either the estimates of the means or the uncertainties, which matches the evolution of the RMS error (Fig. 13). While both $K$ and $\varepsilon$ quickly converge to the values used in the reference simulation (Tab. 4), the low uncertainty in $K$ compared to $\varepsilon$ reflects the physical significance of this parameter: The nonlinear dynamics of a premixed flame are highly sensitive to the timings of pinch-off events [77], which in turn depend on the phase speed at which perturbations travel along with the flame. In our reduced-order model based on the $G$-equation, the phase speed is regulated by the parameter $K$, which makes it highly observable even in the presence of observation noise. Similar to sensitivity analysis [78, 79, Magri2019_amr], this physical insight is gained from the inspection of the uncertainties, which here exceed the residuals between the means and the set of reference values by several orders of magnitude.

Finally, combined state and parameter estimation using the ensemble Kalman smoother is performed (Theorem 6). In Fig. 14, the marginal distributions in the normalized residuals of $K$ and $\varepsilon$ are shown for both the ensemble Kalman filter and smoother. The effect of the ensemble Kalman filter has already been described in Fig. 14. In the forward-backward implementation, the ensemble Kalman smoother takes the last solution of the ensemble Kalman filter and works itself backwards in time (Theorem 2). In terms of information theory, the ensemble Kalman smoother takes at every timestep $k$ all observations, from the past, present and future, into account compared to the ensemble Kalman filter, which only takes the past and the present into account. Statistically speaking, the smoothed distributions are more strictly conditioned than the filtered distributions, over $y_{1:N}$ compared to $y_{1:k}$ (Tab. 3). This surplus in information is evident in the form of lower error bars not just towards the end of the combined state and parameter estimation, but extending all the way to the very first instances of data assimilation. The evaluation of Eq. 6 does not rely on the solution of the governing equations, but instead relies on the predicted and filtered distributions in storage. The ensemble Kalman smoother is thus a computationally inexpensive tool to quantify the uncertainties especially at the beginning of the simulation, where little retrospective data is available. In analogy to direct-adjoint looping [81, 82], combined state and parameter estimation based on filtering and smoothing can be used to obtain otherwise inevitably ad-hoc initial conditions and parameters.

4.3 Effect of observation noise

In a first parameter study, the observation error $\sigma_{\epsilon}$ is varied to investigate the effect of observation noise on the level-set data assimilation framework. In Section 4.2, it has been established that combined state and parameter estimation eventually leads to a statistically stationary state. In Fig. 16, the same RMS error is plotted over time for different observation errors $\sigma_{\epsilon}$. This reveals that the sustained low uncertainties are epistemic in nature [83]: Every reduction in the observation error $\sigma_{\epsilon}$ by one order of magnitude reduces the RMS error by one order of magnitude. It is clear that the observation error poses a lower epistemic bound on how low the uncertainty can be reduced by combined state and parameter estimation.

In Fig. 17, the marginal distributions in the normalized residuals of $K$ and $\varepsilon$ are plotted over time for the different observation errors. The means of the residuals quickly vanish for all three observation errors $\sigma_{\epsilon}$. For the lower observation errors, the error bars increasingly fail to contain the zero residual. Considering that three sigmas correspond to a confidence of 99.7 % under the assumption of normal distributions (Section 2.3), the frequency at which the error bars fail leads to the conclusion that the level-set data assimilation framework underpredicts the uncertainties for low observation errors. This is in agreement with the theoretical analysis of the level-set data assimilation framework for the two-dimensional test case (Section 3.4): Although the individual simulations predict the formation of sharp cusps as the perturbations travel on the respective flame surfaces, the cusps are rounded in the mean of the ensemble (Fig. 6, left). Furthermore, the observable sharpness of a cusp depends on the local number of observations of points on the flame surface. A decrease in the observation error $\sigma_{\epsilon}$ without an increase in the number of observations amounts to a relative decrease in the resolution of the observed cusp (Fig. 8). Therefore, the local values of the eikonal fields of the individual simulations significantly deviate from normal distributions (Fig. 6, right), which leads to incorrect distributions and uncertainties.

4.4 Effect of number of observations

In a second parameter study, subsamples of the observation points are used to investigate the effect of the number of observations on the level-set data asssimilation framework. Furthermore, this serves to illustrate the superiority of our analytic view over the geometric view in data assimilation (Section 3): In the geometric view, the observed interface must be parametric or at least as highly resolved as the grid used in the level-set method to allow the optimal interpolation of a sufficiently large number of points on the predicted interfaces. This restriction does not exist in the analytic view. Hence, the effect of dispersed observation points on the level-set data assimilation framework is also studied here.

In Fig. 18, the RMS error is plotted over time for different randomly sampled subsets of the observation points. In Section 4.2, it has been established that combined state and parameter estimation eventually leads to a statistically stationary state after approximately 10,000 timesteps. With 10 % of the observation points, it takes more than 30,000 timesteps to reach a comparable low-uncertainty plateau. With 2 % of the observation points, the RMS error is still decreasing after 50,000 timesteps while repeatedly failing to correctly predict the pinch-off events for the whole ensemble.

In Fig. 19, the marginal distributions in the normalized residuals of $K$ and $\varepsilon$ are plotted over time for the different subsamples. While the means of the residuals quickly vanish for all three subsamples, the confidence levels for the incomplete subsamples improve less rapidly than for the original sample of observation points. This is in agreement with the theoretical analysis of the level-set data assimilation framework for the two-dimensional test case (Section 3.4): In particular for the 2 %-subsample, certain features in the observed interfaces remain underresolved at the individual timesteps (Fig. 7). Nevertheless, the ensemble at a timestep $k$ representing the filtered probability distribution $p(x_{k},\theta\mid y_{1:k},f)$ (Tab. 3) reflects the knowledge of the observations at all previous timesteps due to the Markov chain property of the probabilistic state space model (Section 2.1). In theory, the lack of resolution in the observations is compensated by more instances of data assimilation to accumulate the same amount of information. The result is a similar statistically stationary state reached at a later time.

5 Conclusion

Reduced-order models based on level-set methods are widely used tools to qualitatively capture and track the nonlinear dynamics of an interface. In this paper, we enhance such a level-set model with a statistical learning technique in order to make the model quantitatively predictive. The statistical learning method is Bayesian, so the uncertainty of the outputs is naturally included in this framework. The main ingredients of the level-set statistical learning are (i) a reduced-order model of an interface based on a level-set method with a law of motion; (ii) the ensemble Kalman filter and smoother; and (iii) external reference data with its uncertainty. The ensemble Kalman filter uses Bayes’ rule as its first principle while assuming normal distributions when data is assimilated. The output of the algorithm is the combined state and parameters estimation, which provides the most likely position of the interface and the model’s parameters given observations. These observations can originate from experimental data or high-fidelity simulations.

To show and verify the algorithm, we perform twin experiments, where the observations are produced by the same code. The assimilation of data from experimental data or high-fidelity simulation requires virtually no modification to the proposed algorithm. The algorithm is applied to three test cases, which are of increasing complexity.

First, a one-dimensional case is presented. We choose the Hamilton-Jacobi generating functions over the whole domain to perform data assimilation, which is called the“analytic approach”. The one-dimensional example shows that the analytic approach overcomes the inaccuracy of the set-theoretic approach, which is based on characteristic functions. Moreover, it is argued that the analytic approach is more straightforward to implement that the geometric approach when moving to higher dimensions.

Second, a two-dimensional, time-independent case is presented, where the position of a sharp corner has to be predicted. It is shown that the level-set data assimilation framework proposed is fully front capturing. The effect that the number of observations has on capturing a sharp corner is investigated. The sensitivity of the analysis solution to each observation is calculated. The one- and two-dimensional cases show that the analytic approach we adopt is more robust than the geometric and set-theoretic approaches that have been used in the past.

Third, the time-dependent nonlinear dynamics of a conical forced premixed flame is studied. The most uncertain states are pinch-offs and the formation of sharp cusps, which are highly nonlinear topological changes of the interface. The two most uncertain parameters are the amplitude of the velocity perturbation at the flame’s base and the perturbation phase speed that wrinkles the flame. In the twin experiment, the combined state and parameter estimation fully recovers the reference solution, which validates the algorithm. Finally, the effect of uncertainty and number of observations is analysed. The uncertainty on the observation most influences the detection of the pinch-off events or the formation of sharp cusps. The data assimilation analytic approach is shown to be versatile because it can assimilate data any time it becomes available, avoiding the cumbersome parameterization required in a geometric approach.

The level-set data assimilation framework is applicable to a number of problems in computational physics involving the motion of an interface [1]. For example, a physical model coupling the dynamics of a premixed flame with an acoustic model of a duct may be used for a data-driven investigation of the Rijke tube [9, 10]. More generally, a number of computational fluid dynamics (CFD) models of premixed combustion, short of direct numerical simulations, combine the laws of motion with conservation laws for the inert flow, e.g. the Bray-Moss-Libby model or the flamelet model [8]. Although it is tempting to directly apply data assimilation to the primitive variables in the Navier-Stokes equations, the discontinuities in density, temperature and species concentrations in light of the theoretical analysis of the set-theoretic approach to level-set data assimilation raise doubts as to whether the kinematics of a premixed flame, which are crucial to thermoacoustics, will be correctly predicted.

Future work will focus on extending the applications of the proposed level-set data assimilation framework. When applied to experiments, the effects of imperfect models will need to be mitigated [66], for example with techniques such as covariance inflation [84, 85] and localization [86]. Future work will also assess the uncertainties in the model, as opposed to uncertainties in the state and parameters.

Acknowledgments

The authors would like to thank the Data Assimilation Research Centre at the University of Reading for the useful discussions at their 2018 Data Assimilation Training. H. Yu is supported the Cambridge Commonwealth, European & International Trust under a Schlumberger Cambridge International Scholarship. L. Magri is supported by the Royal Academy of Engineering under the RAEng Research Fellowships scheme and partially by the visiting fellowship at the Technical University of Munich – Institute for Advanced Study, funded by the German Excellence Initiative and the European Union Seventh Framework Programme under grant agreement n. 291763.

Appendix A Bayesian Inference

Theorem 7** (Marginalization).**

For a joint probability distribution $p(X,Y)$, the marginal probability distributions $p(X)$ and $p(Y)$ are given by

[TABLE]

[TABLE]

Theorem 8** (Bayes’ rule).**

The conditional probability distributions $p(X\mid Y)$ and $p(Y\mid X)$ are related to the joint probability distribution via

[TABLE]

The conditional probability distribution $p(X\mid Y)$ is sometimes referred to as the inverse of the conditional probability distribution $p(Y\mid X)$. The inverse is given by Bayes’ rule:

[TABLE]

If $p(X)$ and $p(Y\mid X)$ are given, $p(Y)$ can be computed through marginalization.

Appendix B Hamiltonian mechanics

Lemma 1** (Lagrangian).**

[TABLE]

Proof.

The generating function is formally the same as the action integral [55, Chapter 9]. Thus, the Lagrangian $\mathcal{L}(r,\dot{r},t)$ is the total derivative of the generating function:

[TABLE]

∎

Lemma 2** (Hamiltonian).**

[TABLE]

Proof.

A Legendre transform of the Lagrangian $\mathcal{L}(r,\dot{r},t)$ (Eq. (76)) gives the Hamiltonian $\mathcal{H}(r,n,t)$:

[TABLE]

It remains to be shown that $\mathcal{H}(r,n,t)$ is indeed a Hamiltonian. The constraints, Eq. (35) and (36), are only functions in $r$ and $n$ respectively. Hence, a suitable choice of Lagrange multipliers $\lambda_{r}$ and $\lambda_{n}$ gives

[TABLE]

[TABLE]

Eq. (80) follows from the laws of motion (Eq. (33)). Eq. (81) follows from tensor calculus:

[TABLE]

where $\mu_{n}$ is the Lagrange multiplier for the constraint imposed by Eq. (36). The gradient of the total derivative of the generating function gives an expression for the partial derivative of the normal vector (Lemma 1):

[TABLE]

The dot product of the normal vector with the total derivative of the normal vector gives an expression for the Lagrange multiplier $\mu_{n}$:

[TABLE]

This gives:

[TABLE]

Note that $\nabla\times n=0$ because $n=\nabla G$. This concludes the proof that $\mathcal{H}(r,n,t)$ is a Hamiltonian. ∎

Theorem 9** (Hamilton-Jacobi equation).**

[TABLE]

Proof.

The Hamilton-Jacobi equation is given by [55, Chapter 9]

[TABLE]

∎

Appendix C Level-set methods

Lemma 3**.**

The covariance matrix of a gradient in a filtered generated function without reinitialization is given by

[TABLE]

where $D$ denotes the gradient operator.

Proof.

Eq. (23) gives

[TABLE]

The probabilistic state space model (Eq. (1), (2)) gives

[TABLE]

[TABLE]

Substituting Eq. 97 into Eq. 96 gives

[TABLE]

∎

Lemma 4**.**

The expected value of a slope in a filtered generated function without reinitialization satisfies

[TABLE]

Proof.

In general, the expected value of a slope is related to the covariance matrix of the corresponding gradient in the following way:

[TABLE]

Substituting Eq. (99) into Eq. (101) gives

[TABLE]

In the first term, the slope $\left|D\psi^{f}\right|$ is equal to 1 due to Huygens’ principle. This gives Lemma 100. ∎

Theorem 10**.**

[TABLE]

Proof.

Applying the following two observations to Eq. (100) gives Theorem 105:

$\mathbb{E}\left[\left|D\psi^{a}\right|\right]^{2}\leq\mathbb{E}\left[\left|D\psi^{a}\right|^{2}\right]$ due to Jensen’s inequality [87]. 2. 2.

$\mathrm{tr}\left(D\left(MC_{\psi\psi}^{f}\right)^{T}\left[C_{\epsilon\epsilon}+MC_{\psi\psi}^{f}M^{T}\right]^{-1}MC_{\psi\psi}^{f}D^{T}\right)\geq 0$ because $D\left(MC_{\psi\psi}^{f}\right)^{T}\left[C_{\epsilon\epsilon}+MC_{\psi\psi}^{f}M^{T}\right]^{-1}MC_{\psi\psi}^{f}D^{T}$ is symmetric and positive semi-definite [88].

∎