2.1. Experimental rig

The experimental rig is shown in Figure 1. It consists of a vertically mounted steel tube with a length of 1 m, inner diameter of 47.4 mm, and wall thickness of 1.7 mm. Both ends of the tube are open to ambient conditions. An electric heating element, shown in detail in Figure 2, is mounted on two support prongs and inserted into the bottom of the tube. The support prongs are attached to an electrically driven traverse so that the heater position can be controlled. The heater is connected to a programmable power supply that is set up to deliver a constant power output. We can therefore study how the thermoacoustic behavior changes with changes in (i) the heater position, and (ii) the heater power.

Figure 1. Diagram of the Rijke tube, rotated for convenience.

Figure 2. (a) Top view, (b) side view, and (c) isometric view of the heater, which consists of two identical concentric annular ceramic plates, each wound with nichrome wire. It is held in place by two threaded support prongs. The dimensions are $ d=47\hskip0.35em \mathrm{mm} $ , $ {d}_i=31.6\hskip0.35em \mathrm{mm} $ , $ {d}_w=0.6\hskip0.35em \mathrm{mm} $ , $ t=5\hskip0.35em \mathrm{mm} $ , $ h=5\hskip0.35em \mathrm{mm} $ , and $ {d}_p=3\hskip0.35em \mathrm{mm} $ . The power is supplied to the nichrome wire by two fabric-insulated copper wires (not shown), which each have diameter $ 4\hskip0.35em \mathrm{mm} $ .

We record the acoustic pressure using eight probe microphones distributed along the length of the tube. The microphones are distributed with 0.1 m spacing, starting at 0.25 m from the bottom of the tube. For a given combination of heater position and heater power, we use a loudspeaker to force the system near its fundamental frequency. We then terminate the forcing and record the pressure time series as the oscillations decay. From this signal, we extract (i) the decay rate, (ii) the natural frequency of oscillations, and (iii) the Fourier-decomposed pressure at each of the microphones, measured relative to a reference microphone. The reference microphone is chosen to be the microphone 0.75 m from the bottom of the tube.

2.2. Thermoacoustic network model

The thermoacoustic oscillations are modeled using a 1D thermoacoustic network model (Dowling and Stow, Reference Dowling and Stow2003; Juniper, Reference Juniper2018). The rig is divided into $ N $ acoustic elements in which forward traveling waves, $ f\left(t-x/c\right) $ , and backward traveling waves, $ g\left(t-x/c\right) $ , propagate. In a given element, the strength of these waves is considered to be constant. The acoustic pressure in the $ i\mathrm{th} $ element is $ {p}_i^{\prime }={f}_i+{g}_i $ , and the acoustic velocity is $ {u}_i^{\prime }=\left({f}_i-{g}_i\right)/\left(\overline{\rho}c\right) $ , where $ \overline{\rho} $ is the local mean density and $ c $ is the local sound speed.

The complex wave amplitudes in adjacent acoustic elements are related through jump conditions for the momentum and energy equations (Juniper, Reference Juniper2018). The jump conditions account for features of the rig that modify the wave amplitudes, such as the visco-thermal dissipation in the boundary layer, the drag and blockage of the heating element, and the fluctuating heat release rate of the heater. The wave reflection at either end of the tube is modeled with complex reflection coefficients. When modeling these jump conditions and boundary conditions, we introduce model parameters, the values of which are not known a priori. As a result, the models tend to be qualitatively accurate, but not quantitatively accurate. We therefore use Bayesian inference to (i) infer the most probable parameter values from data, and (ii) compare several candidate models for the jump conditions, and select the most likely model, given the data.

3.1. Physical system

We consider the task of gathering data by performing experiments on a physical system. The system has a set of design parameters that we control to change the system behavior:

(1) $$ {\mathbf{z}}_i=\mathcal{M}\left({\mathbf{x}}_i\right)+\mathcal{N}\left(0,{\mathbf{C}}_{\mathbf{e}}\right), $$

where $ \mathbf{z} $ is the vector of observed variables, $ \mathcal{M} $ is the measurement operator, and $ \mathbf{x} $ is the vector of design parameters. The observations are corrupted by experimental error, which we assume to be normally distributed with zero mean and covariance $ {\mathbf{C}}_{\mathbf{e}} $ . We therefore only consider the random component of experimental error. Systematic experimental error can be combined with structural error in the model, which can be quantified a posteriori (Yoko and Juniper, Reference Yoko and Juniper2024). We perform a set of $ N $ experiments and use the subscript $ i=\left\{1,2,\dots, N\right\} $ to denote the index of the experiment.

For the Rijke tube described in Section 2.1, the vector of design parameters, $ \mathbf{x} $ , consists of (i) the position of the heater inside the tube, and (ii) the power delivered to the heater. The measurement vector, $ \mathbf{z} $ , contains eight complex numbers: the eigenvalue, whose real part is the decay rate and imaginary part is the natural frequency, and the Fourier-decomposed pressure of seven microphones, measured relative to the reference microphone.

3.2. Physics-based model

We develop several plausible physics-based models that predict the outcome of the experiments:

(2) $$ {\mathbf{s}}_{i,j}={\mathcal{H}}_j\left({\mathbf{x}}_i,\mathbf{a}\right), $$

where $ \mathbf{s} $ are the model predictions, $ \mathcal{H} $ is the candidate physics-based model, $ \mathbf{x} $ is the vector of design parameters, and $ \mathbf{a} $ is a vector of unknown model parameters. We consider $ M $ candidate models and use the subscript $ j=\left\{1,2,\dots, M\right\} $ to denote the index of the candidate model. In the inference framework, we begin by assuming that the model is correct, meaning that we do not explicitly add a term for model error. We expect that, for some set of parameters, $ {\mathbf{a}}_{\mathrm{MP}} $ , the model predictions, $ \mathbf{s} $ , will match the experimental observations, $ \mathbf{z} $ , for all possible design parameters, $ \mathbf{x} $ , within the measurement uncertainty.

For the thermoacoustic network model of the Rijke tube described in Section 2.2, the unknown parameters arise from the sub-models for the visco-thermal dissipation from the boundary layer and the heating element, the reflection coefficients at the ends of the tube, and the fluctuating heat release rate of the heater. We do not know the values of these parameters a priori, and we must often select between several plausible sub-models for each of the physical mechanisms. We address the first problem using Bayesian parameter inference and uncertainty quantification, and the second with Bayesian model comparison.

3.3. Adjoint-accelerated Bayesian inference

3.3.1. Parameter inference

For each candidate model, $ {\mathcal{H}}_j $ , we use the experimental observations to infer the most probable parameters, $ {\mathbf{a}}_{\mathrm{MP}} $ . We begin by assigning a prior probability distribution over the vector of unknown parameters, which we choose to be Gaussian distributed. This prior allows us to incorporate any knowledge we might have about the parameters, while also encoding our confidence in this prior knowledge. We then perform $ N $ experiments to collect the data, $ \mathbf{D}=\left\{{\mathbf{z}}_1,..,{\mathbf{z}}_N\right\} $ . We assimilate the data by performing a Bayesian update on the parameter values:

(3) $$ p\left(\mathbf{a}|\mathbf{D},{\mathcal{H}}_j\right)=\frac{p\left(\mathbf{D}|\mathbf{a},{\mathcal{H}}_j\right)p\left(\mathbf{a}|{\mathcal{H}}_j\right)}{p\left(\mathbf{D}|{\mathcal{H}}_j\right)}. $$

The left-hand side of equation (3), called the posterior, is the probability density of the parameters, given the experimental data and the model. In general, the posterior cannot be evaluated analytically, and numerical computation typically requires millions of model evaluations. This is computationally expensive, even for simple models. At the parameter inference stage, however, we are only interested in finding the most probable parameters, which are those that maximize the posterior. We therefore avoid constructing the full posterior, and instead use an optimization algorithm to find its peak. This is often called the “maximum a posteriori” (MAP) estimate of the parameter values.

The denominator of the right-hand side of equation (3) does not depend on the parameters. We can therefore find the most probable parameters by maximizing the product $ p\left(\mathbf{D}|\mathbf{a},{\mathcal{H}}_j\right)p\left(\mathbf{a}|{\mathcal{H}}_j\right) $ . Recalling that the experimental uncertainty is assumed to be Gaussian distributed, the distribution $ p\left(\mathbf{D}|\mathbf{a},{\mathcal{H}}_j\right) $ is a Gaussian distribution over the data, for a given set of parameters. We also choose the prior distribution, $ p\left(\mathbf{a}|{\mathcal{H}}_j\right) $ , to be Gaussian distributed over the parameters. Under these conditions, we can transform the optimization problem into a quadratic optimization problem by defining the cost function, $ \mathcal{J} $ , to be the negative log of the numerator of equation (3):

(4) $$ \mathcal{J}=\frac{1}{2}\sum \limits_{i=1}^N\left\{{\left({\mathbf{s}}_i\left(\mathbf{a}\right)-{\mathbf{z}}_i\right)}^T{{\mathbf{C}}_{\mathbf{e}}}^{-1}\left({\mathbf{s}}_i\left(\mathbf{a}\right)-{\mathbf{z}}_i\right)+{\left(\mathbf{a}-{\mathbf{a}}_{\mathrm{p}}\right)}^T{{\mathbf{C}}_{\mathbf{a}}}^{-1}\left(\mathbf{a}-{\mathbf{a}}_{\mathrm{p}}\right)+{K}_i\right\}, $$

where $ {\mathbf{a}}_{\mathrm{p}} $ is the vector of prior parameter values, $ {\mathbf{C}}_{\mathbf{a}} $ is the covariance matrix describing the uncertainty in the prior, and $ K $ is a constant from the Gaussian pre-exponential factors, which has no impact on the most probable parameters, $ {\mathbf{a}}_{\mathrm{MP}} $ .

We can solve the optimization problem with the fewest model evaluations by using gradient-based optimization. This requires that we calculate $ \partial \mathcal{J}/\partial \mathbf{a} $ , which we see from equation (4) requires the model parameter sensitivities, $ \partial \mathbf{s}/\partial \mathbf{a} $ . We obtain these using first-order adjoint methods.Footnote ¹ The adjoint method is a technique for calculating the gradient of a function with respect to many parameters, with a computational cost that is independent of the number of parameters (Giles and Pierce, Reference Giles and Pierce2000; Luchini and Bottaro, Reference Luchini and Bottaro2014). The optimization problem can, of course, be solved with any minimization algorithm. However, we will see in the subsequent sections that we will also need the parameter sensitivities for uncertainty quantification, model selection and optimal experiment design. Therefore, it makes sense to use the parameter sensitivities to solve the optimization problem with the fewest model evaluations.

The parameter inference process is illustrated in Figure 3 for a simple system with a single parameter, $ a $ , and a single observable variable, $ z $ . In Figure 3(a), we show the marginal probability distributions of the prior, $ p(a) $ , and the data, $ p(z) $ . The prior and data are independent, so we construct the joint distribution, $ p\left(a,z\right) $ by multiplying the two marginals. In Figure 3(b), we overlay the model predictions, $ s $ , for various values of $ a $ . Marginalizing along the model predictions yields the true posterior, $ p\left(a|z\right) $ . This is possible for a cheap model with a single parameter, but exact marginalization quickly becomes intractable as the number of parameters increases. In Figure 3(c), we plot the cost function, $ \mathcal{J} $ , which is the negative log of the unnormalized posterior. We show the three steps of gradient-based optimization that are required to find the local minimum, which corresponds to the most probable parameters, $ {a}_{\mathrm{MP}} $ .

3.3.2. Uncertainty quantification

Once we have found the most probable parameter values by minimizing equation (4), we estimate the uncertainty in these parameter values using Laplace’s method (Jeffreys, Reference Jeffreys1973; MacKay, Reference MacKay2003; Juniper and Yoko, Reference Juniper and Yoko2022). This method approximates the posterior as a Gaussian distribution with a mean of $ {\mathbf{a}}_{\mathrm{MP}} $ , and a covariance given by the Hessian of the cost function:

(5) $$ {{\mathbf{C}}_{\mathbf{a}}^{\mathrm{MP}}}^{-1}\approx \frac{\partial^2\mathcal{J}}{\partial {a}_l\partial {a}_m}=\sum \limits_{i=1}^N\left\{{{\mathbf{C}}_{\mathbf{a}}}^{-1}+{\mathbf{J}}_i^T{\mathbf{C}}_{\mathbf{e}}^{-\mathbf{1}}{\mathbf{J}}_i+{\left(\mathbf{s}\left(\mathbf{a}\right)-{\mathbf{z}}_i\right)}^T{\mathbf{C}}_{\mathbf{e}}^{-\mathbf{1}}\mathbf{H}\right\}, $$

where the summation is over the $ N $ experimental configurations at which data were collected, which are labeled with the index $ i $ , $ {\mathbf{J}}_i=\mathbf{J}\left({\mathbf{a}}_{\mathrm{MP}},{\mathbf{x}}_i\right) $ is the Jacobian matrix containing the parameter sensitivities of the model predictions, $ \partial {s}_k/\partial {a}_l $ , and $ {\mathbf{H}}_i=\mathbf{H}\left({\mathbf{a}}_{\mathrm{MP}},{\mathbf{x}}_i\right) $ is the rank three tensor containing the second-order sensitivities, $ {\partial}^2{s}_k/\partial {a}_l\partial {a}_m $ . We obtain $ \mathbf{J} $ using first-order adjoint methods, and $ \mathbf{H} $ using second-order adjoint methods.

The accuracy of this approximation depends on the functional dependence between the model predictions and the parameters. This is shown graphically in Figure 4 for three univariate systems. In Figure 4(a), the model is linear in the parameters. Marginalizing a Gaussian joint distribution along any intersecting line produces a Gaussian posterior, so Laplace’s method is exact. In Figure 4(b), the model is weakly nonlinear in the parameters. The true posterior is skewed, but the Gaussian approximation is still reasonable. This panel also shows a geometric interpretation of Laplace’s method: the approximate posterior is given by linearizing the model around $ {\mathbf{a}}_{\mathrm{MP}} $ , and marginalizing the joint distribution along the linearized model. In Figure 4(c), the model is strongly nonlinear in the parameters, so the true posterior is multi-modal and the main peak is highly skewed. In this case, the gradient-based optimization algorithm will only find a single local minimum, which will depend on the choice of initial condition for the optimization. This can be avoided by reducing the extent of the nonlinearity captured by the joint distribution by (i) shrinking the joint distribution by providing more precise prior information or more precise experimental data, or (ii) re-parameterizing the model to reduce the strength of the nonlinearity (MacKay, Reference MacKay2003, Chapter 27).

Figure 3. Illustration of parameter inference on a simple univariate system. (a) The marginal probability distributions of the prior and data, $ p(a) $ and $ p(z) $ , as well as their joint distribution, $ p\left(a,z\right) $ are plotted on axes of parameter value, $ a $ , vs. observation outcome, $ z $ . (b) The model, H, imposes a functional relationship between the parameters, $ a $ , and the predictions, $ s $ . Marginalizing along the model predictions yields the true posterior, $ p\left(a|z\right) $ . This is computationally intractable for even moderately large parameter spaces. (c) Instead of evaluating the full posterior, we use gradient-based optimization to find its peak. This yields the most probable parameters, $ {a}_{\mathrm{MP}} $ .

Figure 4. Illustration of uncertainty quantification for three univariate systems, comparing the true posterior, $ p\left(a|z\right) $ to the approximate posterior from Laplace’s method, $ p{\left(a|z\right)}_L $ . (a) The model is linear in the parameters, so the true posterior is Gaussian and Laplace’s method is exact. (b) The model is weakly nonlinear in the parameters, the true posterior is slightly skewed, but Laplace’s method yields a reasonable approximation. (c) The model is strongly nonlinear in the parameters, the posterior is multi-modal and Laplace’s method underestimates the uncertainty.

In many cases, this approximation will be justifiable, given the substantial reduction in computational cost compared to sampling methods, which are the only viable alternative for constructing the posterior. In previous work, we compared the computational cost of our framework to two sampling approaches (Yoko and Juniper, Reference Yoko and Juniper2024). The comparison was done on a computationally cheap thermoacoustic network model, similar to the one used in this study. Applying our framework to this model, we can compute the posterior probability of five unknown parameters in under 5 s on a single core on a laptop. The same inference problem takes 35 CPU hours running on a workstation when solved with Markov Chain Monte Carlo, and 22 CPU hours when solved with Hamiltonian Monte Carlo.

3.3.3. Uncertainty propagation

We have quantified the uncertainty in the parameter values, but we are also interested in how the parametric uncertainty affects the model predictions. We quantify the prediction uncertainty by propagating the parameter uncertainty through the model. This is done cheaply by linearizing the model around $ {\mathbf{a}}_{\mathrm{MP}} $ and propagating the uncertainties through the linear model:

(6) $$ {\mathbf{C}}_{\mathbf{s}i}={\mathbf{J}}_i^T{\mathbf{C}}_{\mathbf{a}}{\mathbf{J}}_i, $$

where $ {\mathbf{C}}_{\mathbf{s}i} $ is the covariance matrix describing the uncertainty in the model predictions at $ {\mathbf{s}}_i={\mathcal{H}}_j\left({\mathbf{x}}_i,{\mathbf{a}}_{\mathrm{MP}}\right) $ . The marginal uncertainty of each predicted variable is given by the diagonal elements of $ {\mathbf{C}}_{\mathbf{s}i} $ , because the prediction uncertainties are Gaussian.

3.3.4. Model comparison

In many cases, we are faced with a set of candidate models and are required to identify which model is the best. In the Bayesian framework, we can compare candidate models by calculating the posterior probability of each model, given the data:

(7) $$ p\left({\mathcal{H}}_j|\mathbf{D}\right)\hskip0.35em \propto \hskip0.35em p\left(\mathbf{D}|{\mathcal{H}}_j\right)p\left({\mathcal{H}}_j\right). $$

The first factor on the right-hand side of equation (7) is the denominator of equation (3), which is referred to as the marginal likelihood (ML). The second factor is the prior probability that we assign to each model. If we have no reason to prefer one model over another, we assign equal prior probabilities to all models and rank them according to their ML. The ML is calculated by integrating the numerator of equation (3) over parameter space. For a Gaussian posterior, this integral is approximated as

(8) $$ p\left(\mathbf{D}|{\mathcal{H}}_i\right)\approx p\left(\mathbf{D}|{\mathbf{a}}_{\mathrm{MP}},{\mathcal{H}}_i\right)\times p\left({\mathbf{a}}_{\mathrm{MP}}|{\mathcal{H}}_i\right){\left|{\mathbf{C}}_{\mathbf{a}}^{\mathrm{MP}}\right|}^{1/2}. $$

The ML can be broken down into two components, which provide insight into the model comparison. The first factor on the right-hand side of equation (8) is referred to as the best fit likelihood (BFL), which is a measure of how well the model fits the data. The second factor, called the Occam factor (OF), penalizes the model based on its complexity, where the complexity is measured by how precisely the parameter values must be tuned for the model to fit the data to within the experimental uncertainty. A model with a large BFL fits the data well, and vice-versa. A model with a large OF has low parametric complexity, and is unlikely to over-fit the data. Conversely, a model with a small OF is likely to over-fit the data. Therefore, the model with the largest ML is the simplest model that is capable of describing the data, for given measurement error and given priors. The Bayesian model comparison framework therefore naturally enforces Occam’s razor to select the best model.

4.1. Optimal design for parameter inference

We consider the situation where we have collected and assimilated a set of data $ {\mathbf{D}}_i=\left\{{\mathbf{z}}_1,..,{\mathbf{z}}_i\right\} $ , where $ i $ could be zero.Footnote ² We want to know which experiment to perform next in order to gain maximal information about the unknown parameters. In this case, a suitable utility function is the information content of the next experiment, which we can calculate as the change in Shannon entropy between the prior and posterior parameter distributions (Lindley, Reference Lindley1956). The Shannon entropy of a probability distribution is defined as

(10) $$ {S}_i=-{\int}_{\mathbf{A}}P\left(\mathbf{a}|{\mathbf{D}}_i\right){\log}_2\left(p\left(\mathbf{a}|{\mathbf{D}}_i\right)\right)\mathrm{d}\mathbf{a}, $$

where $ {S}_i $ is the Shannon entropy of the parameter probability distribution after the ith experiment has been assimilated.Footnote ³ The utility function, $ u $ , which we have chosen to be the information content of experiment $ i+1 $ , is therefore,

(11) $$ {\displaystyle \begin{array}{c}u\left(\mathbf{x},\mathbf{z},\mathbf{a}\right)=\Delta {S}_{i+1}={S}_i-{S}_{i+1}\\ {}=-{\int}_{\mathbf{A}}P\left(\mathbf{a}|{\mathbf{D}}_i\right){\log}_2\left(p\left(\mathbf{a}|{\mathbf{D}}_i\right)\right)\mathrm{d}\mathbf{a}+{\int}_{\mathbf{A}}P\left(\mathbf{a}|\left\{{\mathbf{D}}_i,{\mathbf{z}}_{i+1}\right\}\right){\log}_2\left(p\left(\mathbf{a}|\left\{{\mathbf{D}}_i,{\mathbf{z}}_{i+1}\right\}\right)\right)\mathrm{d}\mathbf{a}.\end{array}} $$

We note that the utility function involves integration over $ \mathbf{a} $ , so it is no longer a function of the unknown parameters, that is, for this choice of utility function, $ u\left(\mathbf{a},\mathbf{x},\mathbf{z}\right)=u\left(\mathbf{x},\mathbf{z}\right) $ . The expected utility therefore simplifies to

(12) $$ {\displaystyle \begin{array}{c}U={\int}_{\mathbf{Z}}{\int}_{\mathbf{A}}u\left(\mathbf{x},\mathbf{z}\right)p\left(\mathbf{a}|\mathbf{x},\mathbf{z}\right)p\left(\mathbf{z}|\mathbf{x}\right)\mathrm{d}\mathbf{a}\mathrm{d}\mathbf{z}\\ {}={\int}_{\mathbf{Z}}u\left(\mathbf{x},\mathbf{z}\right)p\left(\mathbf{z}|\mathbf{x}\right)\mathrm{d}\mathbf{z}.\end{array}} $$

To evaluate the remaining integral, we once again exploit Laplace’s method, recalling that the probability distribution over the parameters is always Gaussian in our framework. The difference in Shannon entropy between two Gaussians is given by

(13) $$ u\left(\mathbf{x},\mathbf{z}\right)=\Delta {S}_{i+1}=\frac{1}{2}{\log}_2\left(\frac{{\left|{{\mathbf{C}}_{\mathbf{a}}}^{-1}\right|}_{i+1}}{{\left|{{\mathbf{C}}_{\mathbf{a}}}^{-1}\right|}_i}\right) $$

from which we see that maximizing the information gained about the parameters is equivalent to finding the maximum reduction in parameter uncertainty. Recall that, from Laplace’s approximation in equation (5), assimilating the measurement $ {\mathbf{z}}_{i+1} $ updates the parameter uncertainty according to

(14) $$ {\left({{\mathbf{C}}_{\mathbf{a}}}^{-1}\right)}_{i+1}\approx {\left({{\mathbf{C}}_{\mathbf{a}}}^{-1}\right)}_i+{\mathbf{J}}_{i+1}^T{{\mathbf{C}}_{\mathbf{e}}}^{-1}{\mathbf{J}}_{i+1}+{\left(\mathbf{s}{\left(\mathbf{a}\right)}_{i+1}-{\mathbf{z}}_{i+1}\right)}^T{{\mathbf{C}}_{\mathbf{e}}}^{-1}{\mathbf{H}}_{i+1}, $$

where the model predictions, $ \mathbf{s}{\left(\mathbf{a}\right)}_{i+1} $ , Jacobian, $ {\mathbf{J}}_{i+1} $ , and Hessian, $ {\mathbf{H}}_{i+1} $ , are evaluated with the previous most probable parameters, $ {\mathbf{a}}_i $ , and with the candidate experiment design parameters, $ {\mathbf{x}}_{i+1} $ .

From equations (13) and (14), we can see that the choice of which experiment to perform next is dependent on the outcome of that experiment, $ {\mathbf{z}}_{i+1} $ , as one might expect. However, the data-dependence only occurs in the second-order term in equation (14), which is exactly zero for models that are linear in the parameters, and is often small compared to the first-order term for nonlinear models.Footnote ⁴ We therefore proceed with a first-order approximation, in which the posterior covariance is independent of the experimental outcome:

(15) $$ {\left({{\mathbf{C}}_{\mathbf{a}}}^{-1}\right)}_{i+1}\approx {\left({{\mathbf{C}}_{\mathbf{a}}}^{-1}\right)}_i+{\mathbf{J}}_{i+1}^T{{\mathbf{C}}_{\mathbf{e}}}^{-1}{\mathbf{J}}_{i+1}. $$

If the second-order sensitivities are found to be non-negligible for a model of interest, we may still be able to neglect the data-dependent term on the grounds that the discrepancy, $ \left(\mathbf{s}{\left(\mathbf{a}\right)}_{i+1}-{\mathbf{z}}_{i+1}\right) $ , will become small as the model is updated with more data. In this case, we expect that the initial experiments may not be optimal, but that the chosen experiments will become optimal as more experiments are assimilated and the discrepancy decreases. In many cases, this may be an acceptable sacrifice for the computational cost savings of the proposed approach.

By neglecting the data-dependent term, the utility function becomes independent of $ \mathbf{z} $ , and we are able to plan the subsequent experiment using only the model and its adjoint:

(16) $$ U=\Delta {S}_{i+1}=\frac{1}{2}{\log}_2\left(\frac{\left|{\left({{\mathbf{C}}_{\mathbf{a}}}^{-1}\right)}_i+{\mathbf{J}}_{i+1}^T{{\mathbf{C}}_{\mathbf{e}}}^{-1}{\mathbf{J}}_{i+1}\right|}{\left|{\left({{\mathbf{C}}_{\mathbf{a}}}^{-1}\right)}_i\right|}\right). $$

We follow a greedy sequential experiment design process, where we iteratively (i) use equation (16) to identify the most informative experiment, (ii) perform that experiment to collect the data, $ {\mathbf{z}}_{i+1} $ , (iii) minimize equation (4) to find $ {\mathbf{a}}_{i+1} $ , and (iv) solve equation (5) to find $ {\left({\mathbf{C}}_{\mathbf{a}}\right)}_{i+1} $ .

To assist with interpretation of equation (16), we consider the special case of a system with univariate design and observable variables (multiple unknown parameters are still permitted). Under these conditions, the Jacobian, $ \mathbf{J} $ , reduces to a column vector, $ \mathbf{j} $ , and the measurement covariance, $ {\mathbf{C}}_{\mathbf{e}} $ , reduces to the variance, $ {V}_e $ . The quantity $ {\mathbf{J}}^T{{\mathbf{C}}_{\mathbf{e}}}^{-1}\mathbf{J} $ in the numerator of equation (16) reduces to $ {V}_e^{-1}{\mathbf{j}}^T\mathbf{j} $ , which is a rank-one perturbation of the inverse covariance matrix $ {{\mathbf{C}}_{\mathbf{a}}}^{-1} $ . Using the identities for rank-one perturbations, equation (16) reduces to

(17) $$ U=\frac{1}{2}{\log}_2\left(1+{V}_e^{-1}{\mathbf{j}}_{i+1}^T{\mathbf{C}}_{\mathbf{a}i}{\mathbf{j}}_{i+1}\right). $$

This recovers the result that MacKay (Reference MacKay1992) arrived at for scalar interpolation problems: to maximize information gain we must (i) maximize $ {V}_e^{-1} $ , which is to say that we learn the most when we make precise observations, or (ii) maximize $ {\mathbf{j}}^T{\mathbf{C}}_{\mathbf{a}}\mathbf{j} $ , which is the posterior variance of the model predictions, as shown in equation (6). This produces the intuitive result that we learn the most about the unknown parameters when we perform experiments with the design parameters at which our model is most uncertain. This is similar to the optimal sampling approach commonly used in training Gaussian processes, where the acquisition function is chosen to be the variance in the Gaussian process prediction (Huang et al., Reference Huang, Allen, Notz and Zeng2006; Sengupta and Juniper, Reference Sengupta and Juniper2022).

While this aids with interpretation, MacKay’s result is not directly applicable to multivariate systems, unless an ad hoc decision is made about which variable’s uncertainty should drive the decision of which experiment to perform next. This is avoided in our framework by maximizing the change in Shannon entropy, which automatically balances the information gained from each observed variable based on how sensitive it is to the unknown parameters.

4.2. Optimal sensor placement for parameter inference

We now consider the case where one or more of the observed variables, $ \mathbf{z} $ , are spatially varying and are observed with point measurement sensors. We would like to know where to place the sensors in order to gain as much information as possible about the unknown parameters. We may have an existing rig with a fixed number of sensors, and we would like to know where to place the sensors to gain maximal information. Alternatively, we may be designing a new rig, for which we already have a qualitative model, and we would like to know (i) how many sensors we need to buy, and (ii) where we should make provision for instrument access.

To answer these questions, we add the sensor locations to the vector of design parameters, $ \mathbf{x} $ , and find the design parameters that maximize the change in Shannon entropy in the same way as in Section 4.1. We place the sensors sequentially, with each sensor location selected to provide maximum information based on (i) the model and (ii) the information from the existing sensor layout.

This process naturally accounts for the local reduction in information in the vicinity of existing sensors, with the correlation length determined by the model sensitivity. This removes the need to define ad hoc methods to avoid sensor clustering, as done in previous studies on sparse sensor placement (Papadimitriou and Papadimitriou, Reference Papadimitriou and Papadimitriou2015; Bidar et al., Reference Bidar, Anderson, Qin and Bank2024).

4.3. Optimal design for model comparison

Finally, we consider the case where we are trying to identify the best model from a set of candidate models $ {\mathcal{H}}_j,j=\left\{1,2\dots, M\right\} $ . We assume that we have already performed the optimal experiments to learn the unknown parameters of each model. We now want to identify the optimal experiment design, $ {\mathbf{x}}_{i+1} $ , which maximizes the discrimination between the candidate models.

At the (as yet undetermined) experiment design $ {\mathbf{x}}_{i+1} $ , each candidate model will make a slightly different prediction, $ {\mathbf{s}}_{i+1,j} $ , with slightly different uncertainty, $ {\left({\mathbf{C}}_{\mathbf{s}}\right)}_{i+1,j} $ . Before we make the next observation, $ {\mathbf{z}}_{i+1} $ , each model encodes a belief that the data will occur with a probability distribution given by a Gaussian centered around the model prediction:

(18) $$ {\displaystyle \begin{array}{l}{P}_j=p\left({\mathbf{z}}_{i+1}|{\mathcal{H}}_j\right)=\mathcal{N}\left({\mu}_j,{\mathbf{C}}_j\right),\\ {}{\mu}_j={\mathbf{s}}_{i+1,j}={\mathcal{H}}_j\left({\mathbf{x}}_{i+1},{\mathbf{a}}_{i,j}\right),\\ {}{\mathbf{C}}_j={\left({\mathbf{C}}_{\mathbf{s}}\right)}_j+{\mathbf{C}}_{\mathbf{e}},\end{array}} $$

where $ {\mathbf{C}}_{\mathbf{s}} $ is the prediction uncertainty, which is calculated using equation (6), and we have introduced the shorthand $ {P}_j=p\left({\mathbf{z}}_{i+1}|{\mathcal{H}}_j\right) $ to simplify the subsequent notation. This is illustrated for two candidate models in Figure 5.

Figure 5. (a) At each candidate experiment design, $ {x}_{i+1} $ , the two candidate models make slightly different predictions, with different uncertainties. (b) Each model encodes a belief that the next data point will fall within the distribution $ p\left({z}_{i+1}|{H}_j\right) $ .

In order to maximally discriminate between the models, we want to choose $ {\mathbf{x}}_{i+1} $ such that the distributions $ {P}_j $ , where $ j=\left\{1,2,\dots, M\right\} $ , are as dissimilar as possible. In this way, the models that make poor predictions receive the maximum penalty when the data arrive. We therefore choose our utility function to be a measure of dissimilarity between the distributions. A common measure for the dissimilarity between two distributions is the directed Kullback–Leibler (KL) divergence (Kullback and Leibler, Reference Kullback and Leibler1951):

(19) $$ {D}_{KL}\left({P}_j\Big\Vert {P}_k\right)={\int}_{\mathbf{Z}}{P}_j{\log}_2\left(\frac{P_j}{P_k}\right)\mathrm{d}\mathbf{z}, $$

which measures how unexpected the measurement, $ {\mathbf{z}}_{i+1} $ , would be if we used $ {\mathcal{H}}_k $ as our model when $ {\mathcal{H}}_j $ was the better model. This is a suitable measure for our experimental goals, but it only compares two distributions, and we would like to compare multiple models at once. In order to measure the dissimilarity between $ M $ distributions, we use the average divergence (Sgarro, Reference Sgarro1981) as our utility function:

(20) $$ u=\frac{1}{M\left(M-1\right)}\sum \limits_{j=1}^M\sum \limits_{k=1}^M{D}_{KL}\left({P}_j\Big\Vert {P}_k\right), $$

which computes the mean directed KL divergence between all pairs of models. The inner summation measures the average rate at which our confidence in each of the models goes to zero if $ {\mathcal{H}}_j $ is the best model. The outer summation averages these rates while proposing each model as the best model. By maximizing this quantity, we find the experiment design for which the measurement would surprise us the most, thereby reducing our confidence in the largest number of bad models at once.

As before, evaluating the expected utility requires the solution of multiple high-dimensional integrals. Substituting Gaussian distributions for $ {P}_j $ and $ {P}_k $ , and following a similar analysis to Section 4.1, we find that the expected utility reduces to

(21) $$ U=\frac{1}{2M\left(M-1\right)}\sum \limits_{j=1}^M\sum \limits_{k=1}^M\left({\log}_2\left(\frac{\mid {\mathbf{C}}_j\mid }{\mid {\mathbf{C}}_k\mid}\right)-d+ tr\left({\mathbf{C}}_k^{-1}{\mathbf{C}}_j\right)+{\left({\mu}_j-{\mu}_k\right)}^T{\mathbf{C}}_k^{-1}\left({\mu}_j-{\mu}_k\right)\right), $$

where $ d $ is the number of observed variables. We once again find that the integrals, which typically make optimal experiment design expensive, become trivial within our framework. We are therefore able to identify optimal experiments using only the model and its adjoint, and follow a greedy selection process as in the previous sections.

To gain further insight into the expected utility, we consider the case of comparing two candidate models with univariate design and observable variables. Under these conditions, equation (21) reduces to the symmetric KL divergence between a pair of univariate Gaussians:

(22) $$ u=\frac{1}{2}\left\{\left(\frac{1}{V_1}+\frac{1}{V_2}\right){\left({\mu}_1-{\mu}_2\right)}^2+\frac{{\left({V}_1-{V}_2\right)}^2}{V_1{V}_2}\right\}, $$

where $ {\mu}_j $ and $ {V}_j $ are the expected value and variance of the distributions $ p\left({\mathbf{z}}_{i+1}|{\mathcal{H}}_j\right) $ . As in Section 4.1, we recover the result MacKay (Reference MacKay1992) arrived at for scalar interpolation problems: to maximize the discrimination between models we must (i) gather data where the model predictions maximally disagree, measured relative to the confidence in their predictions, and (ii) gather data where the confidence in the models is maximally different. The first result serves to maximize the best-fit-likelihood reward on the model which fits the new data point better. The second result, which is perhaps less intuitive than the first, serves to maximize the Occam penalty on the model with more parametric flexibility.