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

by   Marcos A. Capistrán, et al.

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.



There are no comments yet.


page 1

page 2

page 3

page 4


A constrained pressure-temperature residual (CPTR) method for non-isothermal multiphase flow in porous media

In petroleum reservoir simulation, the standard preconditioner, the Cons...

A block preconditioner for non-isothermal flow in porous media

In petroleum reservoir simulation, the industry standard preconditioner,...

Reconstructing the thermal phonon transmission coefficient at solid interfaces in the phonon transport equation

The ab initio model for heat propagation is the phonon transport equatio...

Nonlinear mixed-dimension model for embedded tubular networks with application to root water uptake

We present a numerical scheme for the solution of nonlinear mixed-dimens...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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 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. , see [6]. 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. [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


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. [6]: 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

  1. for all .

  2. is a linear function in with and reaches the atmospheric pressure value (therefore, is known).

  3. is a right locally analytic function in .

  4. (3)

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


where .

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 [6] 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 [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 denote the coefficients in equation (5) and

its standard deviation, which we assume is the same for all


We model the prior probability density of

using a hierarchical strategy as follows


where we denote . First, we define the indicator function

Next, we model the prior distribution of

using a Gamma distribution with probability density function

The rationale is that both, the Gamma distribution and the hyperparameter

have 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 [2], 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 [4]). 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

Data: Initial point in the parameter space , number of samples
Result: Samples of the posterior distribution (13)
Initialization: Set , ;
for  do
       draw from the probability transition kernel ;
draw from ;
       if  then
             set ;
             set ;
       end if
end for
Algorithm 1 Standard Metropolis Hastings algorithm.

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. [12]

, by means of the unbiased estimator

of the ratio (14). Here and denotes a sample drawn from the conditional probability distribution . Then, equations (10) and (13) allow us to write


Note that the unknown normalizing constants and do not appear in equation (15) and the quotient is well defined if . This gives rise to Algorithm 2

Data: Initial point with , number of samples
Result: Samples of the posterior distribution (13)
Initialization: Set , ;
for  do
       set ;
       draw from the probability transition kernel ;
       if  then
             draw from the conditional distribution ;
draw from ;
             if  then
                   set and ;
             end if
       end if
end for
Algorithm 2 Single Variable Exchange Metropolis Hastings algorithm.

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. 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 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).

Name symbol value and dimension Thermal expansion coefficient Density Specific heat Heat exchange coefficient Inner radius Outer radius Initial temperature Pressure increase rate Final time
Table 1: Parameter setup. We have used parameters and dimensions as described in  [17]

3.1 Sensitivity analysis

Figure 1 illustrates that perturbations in the thermal conductivity coefficient , given by equation (17

), lead to rather small variance in the corresponding temperature

at 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 [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 s.

Figure 1: Uncertainty propagation. Subplots in the top row depict 100 perturbed samples of the thermal conductivity coeffcient, given by equation (17), for and respectively. Subplots in the middle and the bottom rows depict the variance of the numerical simulation of forward mapping (2) acting on the above conductivity coeffcients at and respectively. The smoothing nature of the forward mapping makes it necessary to aquire temperature data at large integration times, e.g. 1000 seconds.

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 [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 ( 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.

(a) Trace plot
(b) Posterior distribution of
(c) True and estimators Middle row depicts 100 samples of the posterior distribution of the thermal conductivity coefficient. At the bottom row are shown the corresponding temperatures evaluated at and respectively. Hierarchical modeling
Figure 2: Example 1. Subplot (a) shows
(a) Trace plot
(b) Posterior distribution of
(c) Truth and estimators
Figure 3: Example 2. Subplot (a) shows

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.


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.


  • [1] J. Adams. Enzyme inactivation during heat processing of food-stuffs. International journal of food science & technology, 26(1):1–20, 1991.
  • [2] J. M. Bardsley. Gaussian markov random field priors for inverse problems. Inverse Problems and Imaging, 7(2):397–416, 2013.
  • [3] J. Cardoso and A. Emery. A new model to describe enzyme inactivation. Biotechnology and bioengineering, 20(9):1471–1477, 1978.
  • [4] 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.
  • [5] 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.
  • [6] 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.
  • [7] 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.
  • [8] 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.
  • [9] 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.
  • [10] V. Isakov. Inverse problems for partial differential equations, volume 127. Springer, 2006.
  • [11] J. P. Kaipio and C. Fox. The bayesian framework for inverse problems in heat transfer. Heat Transfer Engineering, 32(9):718–753, 2011.
  • [12] G. K. Nicholls, C. Fox, and A. M. Watt. Coupled mcmc with a randomized acceptance probability. arXiv preprint arXiv:1205.6857, 2012.
  • [13] 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.
  • [14] 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.
  • [15] 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.
  • [16] 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.
  • [17] 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.
  • [18] 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.
  • [19] J. Wang and N. Zabaras. Hierarchical bayesian models for inverse problems in heat conduction. Inverse Problems, 21(1):183, 2004.
  • [20] A. Zellner. Optimal information processing and bayes’s theorem. The American Statistician, 42(4):278–280, 1988.