Estimating a pressure dependent thermal conductivity coefficient with applications in food technology

TL;DR

This paper presents a novel Bayesian approach using hierarchical Gaussian Markov random fields and a specialized MCMC algorithm to estimate pressure-dependent thermal conductivity in food technology, addressing smoothing effects and sensitivity issues.

Contribution

It introduces a hierarchical GMRF model and a Single Variable Exchange Metropolis-Hastings algorithm for inverse estimation of thermal conductivity under pressure.

Findings

Large integration times improve parameter identification

A signal-to-noise ratio of about 1000 is needed for reliable estimates

The method effectively estimates pressure-dependent thermal conductivity

Abstract

In this paper we introduce a method to estimate a pressure dependent thermal conductivity coefficient arising in a heat diffusion model with applications in food technology. To address the known smoothing effect of the direct problem, we model the uncertainty of the conductivity coefficient as a hierarchical Gaussian Markov random field (GMRF) restricted to uniqueness conditions for the solution of the inverse problem established in Fraguela et al. Furthermore, we propose a Single Variable Exchange Metropolis-Hastings algorithm to sample the corresponding conditional probability distribution of the conductivity coefficient given observations of the temperature. Sensitivity analysis of the direct problem suggests that large integration times are necessary to identify the thermal conductivity coefficient. Numerical evidence indicates that a signal to noise ratio of roughly 1000 suffices to reliably retrieve the thermal conductivity coefficient.

Figures

Tables

Table 1. Table 1: Parameter setup . We have used parameters and dimensions as described in   [ 17 ]

Name	symbol	value and dimension
Thermal expansion coefficient	$\alpha$	$4.217\times {10}^{-4}{K}^{-1}$
Density	$\varrho$	$1000.6\mathrm{kg}{m}^{-3}$
Specific heat	$C_p$	$3780J{\mathrm{kg}}^{-1}{K}^{-1}$
Heat exchange coefficient	$h$	$28W{m}^{-2}{K}^{-1}$
Inner radius	$r_0$	$0m$
Outer radius	$R$	$0.045m$
Initial temperature	$T_0$	$295K$
Pressure increase rate	$\beta$	$120/61\times {10}^6\mathrm{Pa}{s}^{-1}$
Final time	$t_f$	$1000s$

Equations

$$\left\{\begin{array}[]{ll}\varrho C_{p}\dfrac{\partial T}{\partial t}-k(t)\Delta T=\alpha\dfrac{\rm{d}P}{\rm{d}t}(t)T&\mbox{ in }\>B_{R}\times(0,t_{\rm f})\\[8.53581pt] k(t)\dfrac{\partial T}{\partial{\bf n}}=h\left(T^{e}(t)-T\right)&\mbox{ on $\partial B_{R}\times(0,t_{\rm f})$}\\[8.53581pt] T=T_{0}&\mbox{ on }\>B_{R}\times\{0\},\end{array}\right.$$

$$F(k)=T,$$

$${\int }_0^{t_{\mathrm{f}}}k(t)dt\le \varrho C_p\frac{(R-r_0{)}^2}{4}.$$

$$\frac{\mathrm{d}\mathrm{k}}{\mathrm{d}\mathrm{t}}(t)\le \frac{\alpha \beta }{\varrho C_p}\frac{k(t)}{e^{\frac{\alpha \beta }{\varrho C_p}t}-1},t\in [0,t_{\mathrm{f}}].$$

$$k(t)=\exp (u(t)).$$

$$t_j=j\tau ,j=0,1,...,n,n\in \mathbb{N},\tau =\frac{t_{\mathrm{f}}}{n}$$

$$u(t)=\sum _{i=0}^n{u}_i{\phi }_i(t)$$

$$\tau \sum _{i=0}^{n-1}\frac{\exp (u_{i+1})-\exp (u_i)}{u_{i+1}-u_i}\le \varrho C_p\frac{(R-r_0{)}^2}{4}.$$

$$\frac{\mathrm{d}\mathrm{u}}{\mathrm{d}\mathrm{t}}(t)\le f(t),t\in [0,t_{\mathrm{f}}].$$

$$\frac{u_{i+1}-u_{i-1}}{2}\le \frac{\tau }{3}\left(f\left(\frac{t_i+t_{i-1}}{2}\right)+f(t_i)+f\left(\frac{t_{i+1}+t_i}{2}\right)\right),$$

$${\pi }_{\mathbf{\Theta }}(\theta )={\pi }_{\mathbf{\Theta }}(\hat{u},{\sigma }_2)={\pi }_{\hat{\mathbf{U}}|{\mathbf{\Sigma }}_2}(\hat{u}|{\sigma }_2){\pi }_{{\mathbf{\Sigma }}_2}({\sigma }_2){\chi }_Q(\hat{u}),$$

$${\chi }_Q(u)=\{\begin{array}{ll}1& \text{if u(t)=∑i=0nuiφi(t) satisfies eqrefeq:supp1, eqrefeq:supp2 and ui≥u0=log(k(0))},\\ 0& \text{otherwise.}\end{array}$$

$${\pi }_{{\mathbf{\Sigma }}_2}({\sigma }_2)=\frac{1}{\Gamma (a)b^a}{\sigma }_2^{a-1}\exp \left(-\frac{{\sigma }_2}{b}\right).$$

$$\hat{A}=\frac{1}{\tau^{2}}\left(\begin{array}[]{ccccc}2&-1&&&\\ -1&2&-1&&\\ &\ddots&\ddots&\ddots&\\ &&-1&2&-1\\ &&&-1&2\end{array}\right),$$

$$\det (2\pi \hat{\Sigma }{)}^{-1/2}\exp \left(-\frac{1}{2}{\hat{u}}^T\hat{A}\hat{u}\right)$$

$$\hat{\Sigma}=\left(\begin{array}[]{cc}\hat{\Sigma}_{11}&\hat{\Sigma}_{12}\\ \hat{\Sigma}_{21}&\hat{\Sigma}_{22}\\ \end{array}\right),\;\hat{A}=\left(\begin{array}[]{cc}\alpha&v^{T}\\ v&A\\ \end{array}\right).$$

$$\Sigma ={\hat{\Sigma }}_{22}-{\hat{\Sigma }}_{21}{\hat{\Sigma }}_{11}^{-1}{\hat{\Sigma }}_{12}.$$

$$h(u,{\sigma }_2)=\det (2\pi \Sigma {)}^{-1/2}\exp \left(-\frac{1}{2}(u-\mu {)}^{T}A(u-\mu )\right),$$

$${\pi }_{\mathbf{\Theta }}(\theta )=\frac{h(u,{\sigma }_2)}{Z_0({\sigma }_2)}{\pi }_{{\mathbf{\Sigma }}_2}({\sigma }_2){\chi }_Q(u).$$

$$T_{ij}=F(G(u{)}_j{)}_i+{\xi }_{ij}$$

$${\pi }_{T|\Theta }(T|\theta )=\prod _{i,j}\frac{1}{(\sqrt{2\pi }{\sigma }_1{)}^{2n}}\exp \left(-\frac{1}{2{\sigma }_1^2}{\left(T_{ij}-F{\left(G(u{)}_j\right)}_i\right)}^2\right)$$

$${\pi }_{\Theta |\mathbf{T}}(\theta |T)\propto {\pi }_{T|\mathbf{\Theta }}(T|\theta ){\pi }_{\mathbf{\Theta }}(\theta ).$$

$$\alpha =\min \left\{1,\frac{{\pi }_{\Theta |\mathbf{T}}({\theta }^{\prime }|T)q({\theta }^{\prime }|{\theta }_m)}{{\pi }_{\Theta |\mathbf{T}}({\theta }_m|T)q({\theta }_m|{\theta }^{\prime })}\right\};$$

$$E=\frac{{\pi }_{\Theta |\mathbf{T}}({\theta }^{\prime }|T){\pi }_{U|{\mathbf{\Sigma }}_2}(x|{\sigma }_2)q({\theta }^{\prime }|\theta )}{{\pi }_{\Theta |\mathbf{T}}(\theta |T){\pi }_{U|{\mathbf{\Sigma }}_2}(x|{\sigma }_2^{\prime })q(\theta |{\theta }^{\prime })}$$

$$\begin{array}{rl}E& =\frac{{\pi }_{T|\mathbf{\Theta }}(T|{\theta }^{\prime })\frac{h(u^{\prime },{\sigma }_2^{\prime })}{Z_0({\sigma }_2^{\prime })}{\pi }_{{\mathbf{\Sigma }}_2}({\sigma }_2^{\prime }){\chi }_Q(u^{\prime })\frac{h(x,{\sigma }_2)}{Z_0({\sigma }_2)}q({\theta }^{\prime }|\theta )}{{\pi }_{T|\mathbf{\Theta }}(T|\theta )\frac{h(u,{\sigma }_2)}{Z_0({\sigma }_2)}{\pi }_{{\mathbf{\Sigma }}_2}({\sigma }_2){\chi }_Q(u)\frac{h(x,{\sigma }_2^{\prime })}{Z_0({\sigma }_2^{\prime })}q(\theta |{\theta }^{\prime })}\\ & =\frac{{\pi }_{T|\mathbf{\Theta }}(T|{\theta }^{\prime })h(u^{\prime },{\sigma }_2^{\prime }){\pi }_{{\mathbf{\Sigma }}_2}({\sigma }_2^{\prime }){\chi }_Q(u^{\prime })h(x,{\sigma }_2)q({\theta }^{\prime }|\theta )}{{\pi }_{T|\mathbf{\Theta }}(T|\theta )h(u,{\sigma }_2){\pi }_{{\mathbf{\Sigma }}_2}({\sigma }_2){\chi }_Q(u)h(x,{\sigma }_2^{\prime })q(\theta |{\theta }^{\prime })}.\end{array}$$

$$\alpha =\min \left\{1,\frac{{\pi }_{T|\mathbf{\Theta }}(T|{\theta }^{\prime })h(u^{\prime },{\sigma }_2^{\prime }){\pi }_{{\mathbf{\Sigma }}_2}({\sigma }_2^{\prime })h(x,{\sigma }_2^m)q({\theta }^{\prime }|{\theta }_m)}{{\pi }_{T|\mathbf{\Theta }}(T|{\theta }_m)h(u^m,{\sigma }_2^m){\pi }_{{\mathbf{\Sigma }}_2}({\sigma }_2^m)h(x,{\sigma }_2^{\prime })q({\theta }_m|{\theta }^{\prime })}\right\}$$

$$\text{SNR}=\frac{\text{mean}(T)}{{\sigma }_1}=10^3,$$

$$k(t)=\arctan \frac{t}{30}+\mathrm{0.45.}$$

$$k(t)=\frac{1}{2}\left(\arctan \frac{t-450}{150}+3.5\right).$$

Full text

1 Introduction

In this paper we introduce a method to estimate a time dependent thermal conductivity coefficient arising in a heat diffusion model with applications in food technology. Therefore, we are interested in an inverse problem governed by a linear parabolic partial differential equation. According to Infante et al.  [9] “…In high-pressure processes, food is treated with mild temperatures and high pressures in order to inactivate enzymes and preserve as many of its organoleptic properties as possible…”. Consequently, there are several research efforts to model the dynamics of food preservation processes which rely on high pressures away from thermal equilibrium [1, 3, 7, 5, 15, 13, 9, 14, 17, 16]. A related problem is to infer the thermal conductivity coefficient associated to food preservation processes given measurements of the temperature, see [8]. Of particular importance towards the solution of the inverse problem are results on the structural identifiability of the thermal conductivity coefficient, see [6].

On the other hand, inferring a thermal conductivity coefficient from temperature observations is a statistical inference problem. Indeed, it is necessary to model first a forward mapping defined by a well posed problem for a partial differential equation with initial and boundary conditions (e.g. Fraguela et al. [6] give conditions for a mapping from thermal conductivity coefficient to temperature to be well defined). Secondly, we must postulate an observation mapping from state variables to data taking into account a signal to noise ratio. This double tier modeling renders the inverse problem as a statistical inference problem. Related previous work includes boundary heat flux and heat source reconstruction in heat conduction problems using the Bayesian paradigm by Wang and Zabaras [18, 19]. More in general, Kaipio and Fox [11] offer a review of results and challenges in the Bayesian solution of heat transfer inverse problems. We remark that inverse heat conduction problems are difficult given the smoothing nature of the diffusion operator, see Isakov [10], giving rise to a situation where careful modeling and a large signal to noise ratio is necessary in order to be able to solve the inverse problem in a practical setting. Of note, thermocouples, which are data acquisition devices in the setting considered here, have a signal to noise ratio of roughly $10^{3}$ within all the working temperature regime.

In this paper, we rely on the Bayesian paradigm to model the prior information about the parameter to be inferred. The rationale is as follows. In the solution of inverse problems defined by partial differential equations, it is of paramount importance to model prior distributions that are informed with theoretical features of the governing differential equation. On the other hand, Zellner [20] shows that if Bayes theorem is regarded as a learning process in the context of information theory, then the amount of input information is preserved into the output information coded in the posterior distribution. Consequently, in this paper, we construct a prior distribution that incorporates Theorem 14 of [6], giving rise to an unimodal posterior distribution. We analyze the arising posterior distribution using a specially tailored variant of the Metropolis-Hastings algorithm.

The paper is organized as follows: In Section 2 we describe the direct problem that defines the forward mapping, as well as all theoretical results necessary to pose the two tier formulation of the inverse problem as a statistical inference problem. In Section 3 we show our results and discuss our findings. Finally, in Section 4 we reflect upon the reaches and limitations of our approach and offer some perspectives.

2 Problem statement

We shall consider the mathematical modeling of a food preserving method based on high pressure processing away from thermal equilibrium, see Infante et al. [9]. The physical system is partially observed, thus it gives rise to an inverse problem where we want to infer a thermal conductivity coefficient given temperature observations. More precisely, the thermal conductivity coefficient depends on pressure only, e.g. $k=k(P)$, see [6]. The overall goal is identifying $k(P)$ from measurements of the temperature $T$. Of note, data aquisition design renders pressure $P$ as a known and strictly increasing function of time. Therefore, the problem of identifying $k(P)$ is equivalent to identify $k(t)=k(P(t))$.

2.1 Direct problem or forward mapping

According to Smith et al. [17], a mathematical model of the high pressure food preserving method under consideration can be cast as an initial and boundary value problem for a parabolic partial differential equation as follows

[TABLE]

where $B_{R}\subset R^{2}$ is a disk with center $(0,0)$ and radius $R>0$, and $t_{\rm f}>0$ is the final time. Other model components are as follows: $\alpha>0$ is the coefficient of thermal expansion, $\varrho>0$ the density and $C_{p}>0$ the specific heat; $P\in{\cal C}^{1}([0,t_{\rm f}])$ is the pressure at time $t$; $k(t)\geq k_{0}>0$ is the thermal conductivity; $T^{e}$ is the external temperature; $\mathbf{n}$ is the outward unit normal vector at the boundary of $B_{R}$; $h>0$ is a heat exchange coefficient and $T_{0}$ is the initial temperature (assumed to be constant).

Other conditions being equal, problem (1) has a unique (classical) radial solution $T$, and defines a forward mapping

[TABLE]

i.e., given a thermal conductivity coefficient $k=k(t)$, it is possible to evaluate a unique temperature $T$. In the context in inverse problems, we care about conditions on equation (2) to render parameter $k$ uniquely identifiable given data $T$. We shall summarize known uniqueness results for the inverse problem in the next section and use them to pose the inverse problem as a statistical inference problem.

2.2 Inverse problem

The inverse problem that we will consider in the paper is to estimate the thermal conductivity coefficient $k$ given experimental measurements of the temperature $T$. In order to ensure the uniqueness of function $k$, the following context is established in Fraguela et al. [6]: Temperature is measured at two points, one of them on the boundary of the domain and the other one at distance $r_{0}$ of the center point; also, the following hypotheses are assumed to hold

1. H1

$T^{e}(t)\equiv T_{0}$ for all $t\in[0,t_{\rm f}]$. 2. H2

$P$ is a linear function in $[0,t_{\rm f}]$ with $\frac{\rm{d}P}{\rm{d}t}\equiv\beta>0$ and $P(0)$ reaches the atmospheric pressure value (therefore, $k(0)$ is known). 3. H3

$k$ is a right locally analytic function in $[0,t_{\rm f})$. 4. H4

[TABLE] 5. H5

[TABLE]

Note that hypothesis (H5) means an upper estimate of the logarithmic derivative of $k$. Then, taking in account the positivity of $k$, we write

[TABLE]

In order to discretize $u$, we shall use an equidistant grid on interval $[0,t_{\rm{f}}]$

[TABLE]

and the finite element basis of piecewise linear functions $\{\varphi_{i}\}_{i=0}^{n}$ on $[t_{j},t_{j+1}],i=0,1,\ldots,n-1$ such that $\varphi_{i}(t_{j})=\delta_{ij}$ (Kronecker delta). We assume that

[TABLE]

and we want to determine $u_{i}=u(t_{i})$ for $i\geq 1$ ($u_{0}=\log(k(0))$ is known).

From (5), hypothesis (H4) can be discretized as follows

[TABLE]

On the other hand, in order to increase the robustness of our approximations, we weaken restriction (H5) posing it in variational form. First, if we write (3) in terms of $u$ we have

[TABLE]

where $\displaystyle{f(t)=\frac{\alpha\beta}{\varrho C_{p}}\frac{1}{e^{\frac{\alpha\beta}{\varrho C_{p}}t}-1}}$.

Now, multiplying equation (7) by $\varphi_{i}(t),i=1,2,\ldots,n$, integrating by parts and using Simpson’s rule to approximate the definite integrals, hypothesis (H5) becomes

[TABLE]

for $i=1,2,\ldots,n-1$. Theorem 14 in [6] establishes that the inverse problem has an unique solution if the temperature $T$ is solution of problem (1) and hypotheses (H1)-(H5) hold. Let us introduce the following notation. Let us denote by $Q$ the set of functions $u$ such that the corresponding $k(t)$ in equation (4) satisfies $k_{i}\geq k_{0},i=1,2,\ldots,n$ and hypotheses (H1)-(H5) hold, i.e., $u$ satisfies equations (6) and (8). In Subsection 2.3 we shall use the above notation to introduce a two tiered formulation of the inverse problem.

2.3 Statistical inversion

In this Subsection we shall use the optimality property of Bayes formula as a learning process, see Zellner [20], to incorporate Fraguela et al [6] uniqueness theorem into our prior model of the thermal conductivity coefficient, which we discretize as equation (5). Next, we shall correct the prior model with a likelihood function that depends on both, the forward mapping (2) and data. Finally, we propose a variant of the Metropolis-Hastings algorithm to explore the arising posterior distribution.

For the sake of clarity, throughout the current Subsection we shall use bold uppercase letters to denote random variables, while instances of the random variables are denoted as in the direct problem.

Let ${\bf\hat{U}}=({\bf U}_{0},...,{\bf U}_{n})$ denote the coefficients in equation (5) and ${\bf\Sigma}_{2}$ its standard deviation, which we assume is the same for all ${\bf U}_{i}$, $i=0,...,n.$

We model the prior probability density of ${\bf\Theta}=({\bf\hat{U}},{\bf\Sigma}_{2})$ using a hierarchical strategy as follows

[TABLE]

where we denote $\hat{u}=(u_{0},u_{1},\cdots,u_{n})^{T}$. First, we define the indicator function

[TABLE]

Next, we model the prior distribution of $\sigma_{2}$ using a Gamma distribution with probability density function

[TABLE]

The rationale is that both, the Gamma distribution and the hyperparameter $\sigma_{2}$ have support on the positive real numbers, and the Gamma distribution is characterized by two instrumental parameters $a$, $b$, i.e., $\mathbb{E}(\sigma_{2})=ab$, ${\rm var}(\sigma_{2})=ab^{2}$. In the examples below we choose $a=b=1$ such that the prior distribution is monotonically decreasing.

Now, in order to model $\pi_{\bf\hat{U}|\Sigma_{2}}(\hat{u}|\sigma_{2})$, let us consider the following tridiagonal matrix $\hat{A}\in\mathbb{R}^{(n+1)\times(n+1)}$

[TABLE]

where $\tau=\frac{t_{\rm{f}}}{n}$ as before. Of note, $\hat{A}$ is a symmetric positive definite matrix arising if we discretize the negative Laplacian operator with Dirichlet boundary conditions using centered finite differences. Let us denote $\hat{\Sigma}=\hat{A}^{-1}$. Then, according to Bardsley [2], the expression

[TABLE]

defines a Gaussian Markov Random Field (GMRF) with precision matrix $\hat{A}$ (and, covariance matrix $\hat{\Sigma}$). Moreover, we must condition this GMRF to the fact that the first value $u_{0}$ is known. Let us consider $(1,n)$–block partition of matrices $\hat{\Sigma}$ and $\hat{A}$

[TABLE]

Then, our conditional density can be taken as a multivariate normal $N(\mu,\Sigma)$, where the mean $\mu$ is the vector $u_{0}\hat{\Sigma}_{21}\hat{\Sigma}_{11}^{-1}$ and the covariance matrix $\Sigma$ is the Schur complement of matrix $\hat{\Sigma}_{22}$ in $\hat{\Sigma}$, i.e.

[TABLE]

Now, letting $u=(u_{1},...,u_{n})$ we are ready to set $\pi_{\bf\hat{U}|\Sigma_{2}}(\hat{u}|\sigma_{2})=\frac{h(u,\sigma_{2})}{Z(\sigma_{2})}$ where

[TABLE]

and $Z_{0}(\sigma_{2})$ is some unknown normalization constant.

Therefore, our prior probability density in equation (9) can be written as

[TABLE]

Remark. The single variable exchange variant of the Metropolis-Hastings algorithm described below does not require the knowledge of $Z_{0}(\sigma_{2})$ although it depends on $\sigma_{2}$.

Next, we derive a formula to evaluate the likelihood of data $T$ given an instance of the parameter $u=(u_{1},...,u_{n})$, assuming $u_{0}=\log(k(0))$ is known. Let us denote $G(u)_{j}=\exp\left(\sum_{l=0}^{n}u_{l}\varphi_{l}(t_{j})\right)$. We assume that data satisfies

[TABLE]

where data $T_{ij}=T(r_{i},t_{j})$ is a measurement of the temperature at points $r_{i}$ and times $t_{j}$, while $F\left(G(u)_{j}\right)_{i}$ is the forward mapping (2) acting on a instance $u$ of the parameters, while $\xi_{ij}$ is Gaussian noise with mean [math], standard deviation $\sigma_{1}$ and probability distribution $\xi_{ij}\sim\mathcal{N}(0,\sigma_{1}^{2})$. Under the hypothesis of independence of the $\xi_{ij}$, we can approximate the conditional probability distribution $\pi_{{\bf T|\Theta}}(T|\theta)$ of $T$ by the product

[TABLE]

Likelihood (12) is equivalent to $T_{ij}\sim\mathcal{N}\left(F\left(G(u)_{j}\right)_{i},\sigma_{1}^{2}\right)$, where the mean of $T$ is the solution of the direct problem. Given the likelihood (12) and the prior distribution (10), we obtain a model of the posterior distribution of $\Theta$ given $T$ using Bayes identity

[TABLE]

Of note, no analytical expressions are available of the posterior distribution (13), hence we resort to sampling with Markov Chain Monte Carlo Methods. This is achieved by means of the Metropolis-Hastings algorithm.

Let us denote by $q(\theta^{{}^{\prime}}|\theta)$ a proposal for a probability transition kernel (in the numerical examples in Section 3 we shall use the probability transition kernel known as twalk of Christen and Fox [4]). Let us denote $\theta_{m}$ a state in the parameter space and $\theta^{{}^{\prime}}=(u^{{}^{\prime}},\sigma^{{}^{\prime}}_{2})$ a sample drawn from the proposed probability transition kernel. Then, the standard Metropolis-Hastings algorithm is given by Algorithm 1

Of note, quotient (14) has to be calculated through equation (10), which involves the unknown normalizing constant $Z_{0}(\sigma_{2})$. In order to circumvent this problem, we modify the Single Variable Exchange algorithm described in Nicholls et al. [12], by means of the unbiased estimator

[TABLE]

of the ratio (14). Here $\theta=(u,\sigma_{2})$ and $x$ denotes a sample drawn from the conditional probability distribution $\pi_{{\bf U|\Sigma_{2}}}(\cdot|\sigma_{2}^{\prime})$. Then, equations (10) and (13) allow us to write

[TABLE]

Note that the unknown normalizing constants $Z_{0}(\sigma_{2})$ and $Z_{0}(\sigma_{2}^{\prime})$ do not appear in equation (15) and the quotient is well defined if $u\in Q$. This gives rise to Algorithm 2

Algorithm 2 belongs to a wide class of Metropolis-Hastings algorithms with randomized acceptance probability discussed by Nicholls et al. [12], where the ratio of two target distributions is not available, but the ratio can be approximated by an unbiased estimator, e.g. $E$ in equation (15). The equilibrium distribution of the modified algorithm is the original target distribution.

3 Numerical experiments

In the first part of this Section we explore numerically how the forward mapping propagates uncertainty in the thermal conductivity coefficent. Later, we explore the posterior distribution (13) using Algorithm 2 in two numerical examples. Hypotheses (H1)-(H5) are assumed to hold in all cases. We have used Tylose parameters, a setup of parameter values and dimensions described in Table 1. This parameter set has already been used as a valid synthetic setting in [17]. We have solved the direct problem (1) approximating the polar Laplacian and the derivatives in the boundary conditions by means of standard centered order two finite differences and using the classical Crank-Nicolson time stepping, resulting in a second order method, both in time and space. We have discretized the spatial domain with 102 grid points, and we have taken 350 timesteps. Synthetic data is created solving the same direct problem with a grid 100 times more fine in time. Observations of temperature at both, the center and a point of the boundary of the ball, are assumed to be poluted with noise according to likelihood (12), where we have chosen the noise standard deviation $\sigma_{1}$ such that the signal to noise ratio (SNR) satisfies

[TABLE]

and $\text{mean}(T)$ is the mean of the temperature $T=T(r,t)$ for $(t,r)\in B_{R}\times(0,t_{\rm f})$ measured in Kelvin degrees. We remark that thermocouples used in real applications to aquire temperature data have roughly the same signal to noise ratio (16).

3.1 Sensitivity analysis

Figure 1 illustrates that perturbations in the thermal conductivity coefficient $k=k(t)$, given by equation (17), lead to rather small variance in the corresponding temperature $T=T(r,t)$ at $r=0$ and $r=R$. Left column has signal to noise ratio $SNR=10$ and right column has signal to noise ration $SNR=10^{3}$. This feature of the forward mapping (2) pinpoints how difficult it is to solve the corresponding inverse problem. Indeed, as Kaipio and Fox [11] establish, such a narrow posterior distribution might give rise to a situation where a computational forward mapping leads to an unfeasible posterior model, e.g., a posterior model such that the true thermal conductivity coefficient has arbitrarily small probability. Of note, Figure 1 also illustrates that the temperature variance is an increasing function of time. Therefore it is possible to design observation times large enough to circumvent the smoothing effect of the forward problem. In the examples shown below we have chosen final time $t_{f}=1000$ s.

3.2 Parameter estimation

In each case we have carried out $5\times 10^{5}$ steps of the Markov Chain Monte Carlo. In Algorithm 2, we have used the proposal $q(\theta^{{}^{\prime}}|\theta)$ given by the twalk of Christen and Fox [4]. Our choice was based on the facts that the twalk proposal has no tunning parameters and commutes with affine transformations of space. All programs for the numerical experiments below were coded in python 3.7 and are available in github (https://github.com/MarcosACapistran/diffusion) Following standard notation, we denote by $\theta_{MAP}$ the maximum a posteriori probability estimate and by $\theta_{CM}$ the conditional mean estimate.

Example 1

In the first synthetic example we have used parameters described in Table 1 and thermal conductivity is given by the function

[TABLE]

Example 2

In the second synthetic example, we have kept the same parameters from the previous example, except the thermal conductivity, taking in this case

[TABLE]

Figures 2 and 3 show the results of identyfing the thermal conductivity coefficient $k(t)$ for examples 1 and 2 respectively. Figures 2a and 3a are numerical evidence of the convergence of the Markov Chain Monte Carlo in both examples. We have taken $SNR=10^{3}$. Of note, Figures 2b and 3b show roughly the same marginal posterior distribution for the variance of the prior model for both examples. Finally, Figures 2c and 3c provide further evidence that of the smoothing property of the numerical simulation of the forward mapping (2). Although we have used a low dimensional representation (4)-(5) with $n=10$ in both cases, it is apparent that the true value of the of the thermal conductivity coeffcient lies in the support of the posterior model.

4 Conclusions

In this paper we have shown a method to overcome the typical narrowness of the posterior distribution arising in an inverse heat diffusion problem with applications in food technology. Our strategy is to construct a hierarchical prior model of the thermal conductivity coefficient restricted to uniqueness conditions for the solution of the inverse problem. An important feature of our apprach is that the variance of the prior model is a parameter to be inferred from data. Finally, we propose a Single Variable Exchange Metropolis-Hastings algorithm to correctly sample the arising posterior distribution. Numerical evidence indicates that the resulting posterior model of the thermal conductivity coeffcient contains the true value.

In the numerical implementation, we have resorted to approximation methods to turn the inference of the thermal conductivity coefficient into a parametric problem. We have used a piecewise linear function to approximate the quantity of interest, and inference is carried out over a finite number of real coefficients.

In the narrative of food technology, the results obtained in this paper might serve as a basis to explore the dimension of the data informed subspace as a function of the signal to noise ratio. The reliable solution of the inverse problem with quantifiable uncertainty provides a basis for industrial analysis towards better controled food preservation methods.

Acknowledgements

M. A. C. acknowledges partially founding from ONRG, RDECOM and CONACYT CB-2017-284451 grants. J. A. I. R. acknowledges the financial support of the Spanish Ministry of Economy and Competitiveness under the project MTM2015-64865-P and the research group MOMAT (Ref.910480) supported by Universidad Complutense de Madrid.