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. “…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 . Of particular importance towards the solution of the inverse problem are results on the structural identifiability of the thermal conductivity coefficient, see .
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.  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  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 , 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 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 
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, 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. . 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. , see . The overall goal is identifying from measurements of the temperature . Of note, data aquisition design renders pressure as a known and strictly increasing function of time. Therefore, the problem of identifying is equivalent to identify .
2.1 Direct problem or forward mapping
According to Smith et al. , 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
where is a disk with center and radius , and is the final time. Other model components are as follows: is the coefficient of thermal expansion, the density and the specific heat; is the pressure at time ; is the thermal conductivity; is the external temperature;
is the outward unit normal vector at the boundary of; is a heat exchange coefficient and is the initial temperature (assumed to be constant).
Other conditions being equal, problem (1) has a unique (classical) radial solution , and defines a forward mapping
i.e., given a thermal conductivity coefficient , it is possible to evaluate a unique temperature . In the context in inverse problems, we care about conditions on equation (2) to render parameter uniquely identifiable given data . 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 given experimental measurements of the temperature . In order to ensure the uniqueness of function , the following context is established in Fraguela et al. : Temperature is measured at two points, one of them on the boundary of the domain and the other one at distance of the center point; also, the following hypotheses are assumed to hold
for all .
is a linear function in with and reaches the atmospheric pressure value (therefore, is known).
is a right locally analytic function in .
Note that hypothesis (H5) means an upper estimate of the logarithmic derivative of . Then, taking in account the positivity of , we write
In order to discretize , we shall use an equidistant grid on interval
and the finite element basis of piecewise linear functions on such that (Kronecker delta). We assume that
and we want to determine for ( is known).
From (5), hypothesis (H4) can be discretized as follows
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 we have
Now, multiplying equation (7) by , integrating by parts and using Simpson’s rule to approximate the definite integrals, hypothesis (H5) becomes
for . Theorem 14 in  establishes that the inverse problem has an unique solution if the temperature is solution of problem (1) and hypotheses (H1)-(H5) hold. Let us introduce the following notation. Let us denote by the set of functions such that the corresponding in equation (4) satisfies and hypotheses (H1)-(H5) hold, i.e., 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 , to incorporate Fraguela et al  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 denote the coefficients in equation (5) and
its standard deviation, which we assume is the same for all,
We model the prior probability density ofusing a hierarchical strategy as follows
where we denote . First, we define the indicator function
Next, we model the prior distribution of
The rationale is that both, the Gamma distribution and the hyperparameterhave support on the positive real numbers, and the Gamma distribution is characterized by two instrumental parameters , , i.e., , . In the examples below we choose such that the prior distribution is monotonically decreasing.
Now, in order to model , let us consider the following tridiagonal matrix
where as before. Of note, 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 . Then, according to Bardsley , the expression
defines a Gaussian Markov Random Field (GMRF) with precision matrix (and, covariance matrix ). Moreover, we must condition this GMRF to the fact that the first value is known. Let us consider –block partition of matrices and
Then, our conditional density can be taken as a multivariate normal , where the mean is the vector and the covariance matrix is the Schur complement of matrix in , i.e.
Now, letting we are ready to set where
and is some unknown normalization constant.
Therefore, our prior probability density in equation (9) can be written as
Remark. The single variable exchange variant of the Metropolis-Hastings algorithm described below does not require the knowledge of although it depends on .
Next, we derive a formula to evaluate the likelihood of data given an instance of the parameter , assuming is known. Let us denote . We assume that data satisfies
where data is a measurement of the temperature at points and times , while is the forward mapping (2) acting on a instance of the parameters, while is Gaussian noise with mean , standard deviation and probability distribution . Under the hypothesis of independence of the , we can approximate the conditional probability distribution of by the product
Likelihood (12) is equivalent to , where the mean of 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 given using Bayes identity
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 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 ). Let us denote a state in the parameter space and 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 . In order to circumvent this problem, we modify the Single Variable Exchange algorithm described in Nicholls et al. 
, by means of the unbiased estimator
Algorithm 2 belongs to a wide class of Metropolis-Hastings algorithms with randomized acceptance probability discussed by Nicholls et al. , where the ratio of two target distributions is not available, but the ratio can be approximated by an unbiased estimator, e.g. 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 . 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 such that the signal to noise ratio (SNR) satisfies
and is the mean of the temperature for 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
), lead to rather small variance in the corresponding temperatureat and . Left column has signal to noise ratio and right column has signal to noise ration . This feature of the forward mapping (2) pinpoints how difficult it is to solve the corresponding inverse problem. Indeed, as Kaipio and Fox  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 s.
3.2 Parameter estimation
In each case we have carried out steps of the Markov Chain Monte Carlo. In Algorithm 2, we have used the proposal given by the twalk of Christen and Fox . 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 the maximum a posteriori probability estimate and by 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
- 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
Figures 2 and 3 show the results of identyfing the thermal conductivity coefficient for examples 1 and 2 respectively. Figures (a)a and (a)a are numerical evidence of the convergence of the Markov Chain Monte Carlo in both examples. We have taken . Of note, Figures (b)b and (b)b show roughly the same marginal posterior distribution for the variance of the prior model for both examples. Finally, Figures (c)c and (c)c 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 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.
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.
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.
-  J. Adams. Enzyme inactivation during heat processing of food-stuffs. International journal of food science & technology, 26(1):1–20, 1991.
-  J. M. Bardsley. Gaussian markov random field priors for inverse problems. Inverse Problems and Imaging, 7(2):397–416, 2013.
-  J. Cardoso and A. Emery. A new model to describe enzyme inactivation. Biotechnology and bioengineering, 20(9):1471–1477, 1978.
-  J. A. Christen, C. Fox, et al. A general purpose sampling algorithm for continuous distributions (the t-walk). Bayesian Analysis, 5(2):263–281, 2010.
-  S. Denys, A. M. Van Loey, and M. E. Hendrickx. A modeling approach for evaluating process uniformity during batch high hydrostatic pressure processing: combination of a numerical heat transfer model and enzyme inactivation kinetics. Innovative Food Science & Emerging Technologies, 1(1):5–19, 2000.
-  A. Fraguela, J. Infante, A. Ramos, and J. Rey. A uniqueness result for the identification of a time-dependent diffusion coefficient. Inverse Problems, 29(12):125009, 2013.
-  M. Hendrickx, L. Ludikhuyze, I. Van den Broeck, and C. Weemaes. Effects of high pressure on enzymes related to food quality. Trends in Food Science & Technology, 9(5):197–203, 1998.
-  J. Infante, M. Molina-Rodríguez, and Á. Ramos. On the identification of a thermal expansion coefficient. Inverse Problems in Science and Engineering, 23(8):1405–1424, 2015.
-  J. A. Infante, B. Ivorra, A. Ramos, and J. M. Rey. On the modelling and simulation of high pressure processes and inactivation of enzymes in food engineering. Mathematical Models and Methods in Applied Sciences, 19(12):2203–2229, 2009.
-  V. Isakov. Inverse problems for partial differential equations, volume 127. Springer, 2006.
-  J. P. Kaipio and C. Fox. The bayesian framework for inverse problems in heat transfer. Heat Transfer Engineering, 32(9):718–753, 2011.
-  G. K. Nicholls, C. Fox, and A. M. Watt. Coupled mcmc with a randomized acceptance probability. arXiv preprint arXiv:1205.6857, 2012.
-  T. Norton and D.-W. Sun. Recent advances in the use of high pressure as an effective processing technique in the food industry. Food and Bioprocess Technology, 1(1):2–34, 2008.
-  L. Otero, B. Guignon, C. Aparicio, and P. Sanz. Modeling thermophysical properties of food under high pressure. Critical reviews in food science and nutrition, 50(4):344–368, 2010.
-  L. Otero, A. Ramos, C. De Elvira, and P. Sanz. A model to design high-pressure processes towards an uniform temperature distribution. Journal of Food Engineering, 78(4):1463–1470, 2007.
-  V. Serment-Moreno, G. Barbosa-Cánovas, J. A. Torres, and J. Welti-Chanes. High-pressure processing: kinetic models for microbial and enzyme inactivation. Food Engineering Reviews, 6(3):56–88, 2014.
-  N. Smith, S. Mitchell, and A. Ramos. Analysis and simplification of a mathematical model for high-pressure food processes. Applied Mathematics and Computation, 226:20 – 37, 2014.
J. Wang and N. Zabaras.
A bayesian inference approach to the inverse heat conduction problem.International Journal of Heat and Mass Transfer, 47(17-18):3927–3941, 2004.
-  J. Wang and N. Zabaras. Hierarchical bayesian models for inverse problems in heat conduction. Inverse Problems, 21(1):183, 2004.
-  A. Zellner. Optimal information processing and bayes’s theorem. The American Statistician, 42(4):278–280, 1988.