An inverse estimation of surface temperature using the maximum entropy method

International Communications in Heat and Mass Transfer 34 (2007) 37 – 44 www.elsevier.com/locate/ichmt

An inverse estimation of surface temperature using the maximum entropy method☆ Sun Kyoung Kim a,⁎, Bup Sung Jung b , Woo Il Lee b b

a Department of Die and Mold Design Seoul National University of Technology Seoul, 139-743 Republic of Korea Department of Mechanical and Aerospace Engineering, Seoul National University, Seoul, 151-742, Republic of Korea

Available online 27 September 2006

Abstract The maximum entropy method (MEM) is applied to estimation of surface temperature from temperature readings. The inverse heat conduction problem is reformulated for MEM and a three-phase solution method utilizing the successive quadratic programming (SQP) is addressed. Computational results by the proposed MEM are presented and compared with results by the conventional methods. © 2006 Elsevier Ltd. All rights reserved. Keywords: Inverse heat conduction problem (IHCP); Maximum entropy method (MEM); Adjoint problem; Surface temperature

1. Introduction Retrieving imperfectly described boundary conditions in heat conduction problems from the temperature measurement inside or on the boundary can be facilitated by solving the corresponding inverse heat conduction problem (IHCP) [1]. In general, solving an IHCP is not an easy task since it is mathematically ill-posed. In order to overcome such difficulty and to obtain a reliable inverse solution, a variety of numerical techniques have been proposed. These techniques include the least-squares method with regularization [2], the sequential function specification method [1], space marching techniques [3], and the gradient method utilizing adjoint formalism [4]. A number of modifications to these methods have been made to enhance the solution accuracy and stability. Despite such efforts, conventional methods mentioned above are limited in improving the resolution [5]. In this work, the maximum entropy method (MEM) is employed to estimate the unknown temperature accurately from the measured temperature. The MEM, which is based on probabilistic theory, allows seeking the most likely inverse solution [6]. The MEM has been applied to various ill-posed problems including image reconstruction, ill-posed partial differential equation and economic data recovery. Generally, an inverse problem can be treated as an optimization problem. In the MEM, the consistency between the measurement and the estimation of temperature acts as a constraint of the optimization problem. The objective function of this optimization problem is the information entropy that is defined by Shannon [7]. Jaynes proposed to use the information entropy in the field of statistical mechanics [8]. The ☆

Communicated by: K. Suzuki and S. Nishio. ⁎ Corresponding author. E-mail addresses: [email protected] (S.K. Kim), [email protected] (B.S. Jung), [email protected] (W.I. Lee).

0735-1933/$ - see front matter © 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.icheatmasstransfer.2006.08.011

38

S.K. Kim et al. / International Communications in Heat and Mass Transfer 34 (2007) 37–44

Nomenclature f f F L N S t T x Y Z

Unknown temperature Vector form of f Total sum of the unknown temperature Lagrangian function Number of measurement data Square norm of the temperature difference Time Temperature Distance from exposed surface Temperature reading Conjugate direction vector

Symbols β Step length during CGM γ Conjugate coefficient Δ Variation due to variation of f to f + Δf Δf Variation of unknown temperature ε Iteration stoppage criterion for small σ λ Adjoint variable Λ Lagrange multiplier σ Standard deviation of Y τ ≡ tf − t ∇f Frechet derivative with respect to f Subscripts f Final t ≡ ∂ / ∂t x ≡ ∂ / ∂x xx ≡ ∂2 / ∂x2 Superscript p Denotes the estimated value of the previous iteration

MEM allows obtaining maximum possible information from given measurement data of limited accuracy. That is, if an inverse solution is a maximum entropy (ME) solution, it is the most likely solution among the candidates of many inverse solutions which are consistent with the measurement data. The entropy function has a logarithmic form whereas the constraint has a quadratic form. Consequently, the solution procedure requires a nonlinear constrained optimization, which usually consumes considerable computational resources. In the present study, the successive quadratic programming (SQP) is employed as a solution method for the optimization problem. The SQP is a general and straightforward method applicable for nonlinear constrained optimization [9]. Not many works about IHCP methods based on ME principle and it is still a maturing method. Ramos and Giovannini applied ME principle to detect defects by estimating thermal conductance using tomographic technique [10]. Lair et al. estimated heat flux in die forging using MEM [11]. Muniz et al. solved a retrospective IHCP [4], whose unknown is an initial condition, using MEM and Tikhonov regularization method and they compared the results by both methods [12]. The ME solutions in these studies do not show noticeable improvement. The aforementioned three works do not reveal any characteristic features of MEM such as super-resolution and entropy concentration. Moreover those works do not address the procedure in determining the total sum of unknown components that should be provided as an observed value for ME estimation [13]. Performing simply bounded unconstrained optimization with a fixed

S.K. Kim et al. / International Communications in Heat and Mass Transfer 34 (2007) 37–44

39

Fig. 1. Schematic of considered IHCP.

regularization parameter makes them fail to obtain a true ME solution. That is, by such a method, maximum state of the information entropy cannot be attained. Recently, Kim and Lee proposed a ME-based solution method for onedimensional IHCP, which seeks unknown heat flux history [14]. In the referenced work, several generic features of MEM including remarkable resolution improvement were found. In this work, a typical IHCP in one-dimensional slab with a finite thickness is considered. The unknown temperature of interest is estimated on one boundary. At the same time, temperature readings are acquired on the other boundary which is thermally insulated. The maximum entropy estimation of the time-varying surface temperature is conducted for varying measurement error levels and sampling rates to verify the validity of the developed method. Furthermore, this work proposes a new total sum estimation technique for MEM to improve the accuracy. The computational results of the proposed method based on the maximum entropy principle are compared with the conventional methods in terms of the accuracy of the solution for several stringent cases. 2. Formulation Here, an inverse method for estimating unknown temperature from temperature reading is formulated. MEM demands IHCP to be reformulated into a nonlinear constrained optimization problem, which forms the main structure of the current approach. The optimization is accomplished by a procedure composed of three phases, which are feasibility-achieving phase, quadratic regularization (QR) phase, and maximum entropy phase. In each phase, the objective function, constraints, and their derivatives should be evaluated. The most vital part of this evaluation is conducted basically by solving the partial differential equations such as the direct problem and the adjoint problem. The whole mathematical expressions involved in this method are given below specifically. Temperature, time, and position variables are dimensionless in the following formulation. The domain of interest is a finite one-dimensional slab with length 1 as illustrated in Fig. 1. The temperature of the left side f (t) is unknown and the right side is insulated. Material properties are considered constant over the entire domain. The measured temperature Y(t) is acquired on the insulated surface, x = 1. The governing equation of this transient conduction problem with unit thermal diffusivity is given as follows. Tt ¼ Txx ; T ð0; tÞ ¼ f ðtÞ; Tx ð1; tÞ ¼ 0; T ðx; 0Þ ¼ 0

ð1Þ

When the unknown temperature f is perturbed by a small amount Δf, T undergoes a perturbation T + ΔT. Then, substituting in Eq. (1) T by T = ΔT, f by f + Δf, gives the following sensitivity equation. ðDT Þt ¼ ðDT Þxx ; DT ð0; tÞ ¼ Df ðtÞ; DTx ð1; tÞ ¼ 0; DT ðx; 0Þ ¼ 0

ð2Þ

The current inverse problem aims to achieve an optimal consistency between the measured and computed temperature. Such consistency can be evaluated by the following residual norm. Z S¼ 0

tf

½T ð1; t; f Þ−Y ðtÞ2

ð3Þ

40

S.K. Kim et al. / International Communications in Heat and Mass Transfer 34 (2007) 37–44

In order to derive the adjoint problem, the above expression is modified by involving Eq. (1) as a constraint and introducing a Lagrange multiplier λ as follows. Z tf Z tf Z 1 S¼ ½T ð1; t; f Þ−Y ðtÞ2 dt þ k½Txx −Tt dxdt ð4Þ 0

0

0

Considering ΔS = S( f + Δf ) − S( f R) and neglecting the terms with the second order of smallness, we obtain Rt t R1 DS ¼ 0f 2½T ð1; tÞ−Y ðtÞDð1; tÞdt þ 0f 0 k½ðDT Þxx −ðDT Þx dxdt. Integrating by parts and utilizing the arbitrariness of the increments of state variable ΔT and the boundary and initial conditions of the sensitivity equation, we obtain the following adjoint problem from the above expression −kt ¼ kxx ; kð0; tÞ ¼ 0; kx ð1; tÞ ¼ 2½T ð1; t; f Þ−Y ðtÞ; kðx; tf Þ ¼ 0

ð5Þ

TheR above equation can be solved without difficulty by replacing t by τ where τ = tf − t. Then, ΔS becomes t DS ¼ 0f kx ð0; tÞDT ð0; tÞdt. On the other hand, Rin the square-integrable function space (L2 functional space), the t following expression is valid by definition DS ¼ 0f jf SDT ð0; tÞdt. From the comparison of these two equations, the following gradient equation is obtained. jf S ¼ kx ð0; tÞ

ð6Þ

Generally, regularization is central to an inverse algorithm. A three-phase strategy is employed here to efficiently conduct suitable regularization without resolution loss. The first phase performs the iterative regularization typically done in the conventional gradient method [4]. The second phase performs the conventional regularization (first order, quadratic) [2]. Then, the third phase enhances the accuracy of the solution by the maximum entropy method (MEM). It should be emphasized that the maximum entropy estimation cannot be performed without a prerequisite value that is an integrated value of the unknown over the entire time domain, which is given by Z F¼

tf

f ðtÞdt

ð7Þ

0

This value is usually not available as a priori information in IHCP. Thus, the unknown temperature f(t) for the evaluation of F is obtained by the conventional gradient method. Thanks to the principle of energy conservation, F can be estimated accurately despite the deterministic bias inherent in the inverse estimator [15]. Therefore, the current formalism inevitably comprises multiple phases. The first phase is unconstrained optimization whereas the second and third phases are constrained optimization. The first phase aids the second phase to start with satisfying the feasibility condition of the optimization problem [9]. The conjugate gradient method is adopted for the first phase since it can achieve the feasibility condition efficiently. The second phase fulfills the optimal regularization to calculate F exactly and to provide a good initial guess for the third phase. The conventional gradient method achieves desired regularization by stopping the iteration timely with the aid of its viscous nature [4]. It aims to satisfy the following equality. L1 ð f Þ ¼ S−tf r2 ¼ 0

ð8Þ

Let us call an inverse solution that fulfills the above condition as a feasible solution. In this study, such solution is achieved in the first phase by the conjugate gradient method (CGM) combined with the bisection search. That is, after satisfying an inequality S b tf σ2 using CGM, then a feasible solution is sought by the bisection search, which starts with the initial interval given by the last two solutions of the conjugate gradient search. The violation level of Eq. (8) pffiffiffiffi is set equal to er2 tf = N where e ¼ 0:001 and N is the number of discrete temperature readings. Besides, when the number of iteration exceeds its maximum, which is set as 4N, with failing to achieve the feasibility, the estimation continues in the next phase. The CGM for IHCP is well described for various classes of problems in a previous work [16]. However, it is worth repeating in order to emphasize the specific features of CGM for IHCP. The CGM updates the new unknown by f = f p − βS where p denotes the estimated value of the previous iteration, and β is the step length that minimizes S

S.K. Kim et al. / International Communications in Heat and Mass Transfer 34 (2007) 37–44

41

( f p − βZ ), and p is given by Z = ∇f S + γZ p. Here, γ is the conjugate coefficient that is R tf the conjugate R tf p direction 2 2 calculated R t by g ¼ 0 f ðtÞ dt = 0 f ðtÞ dt. The R t step length can be evaluated using the first order Taylor series expansion as b ¼ 0l ½T ð1; t; f Þ−Y ðtÞDT ð1; t; Df Þdt = 0f DT ð1; t; Df Þ2 dt. Here, ΔT(1, t; Δf ) is obtained by solving the sensitivity equation with setting Δf equal to Z. Then, utilizing the feasible solution from the first phase as an initial guess, the second phase performs the quadratic regularization (QR) by directly achieving the desired state of the inverse solution, which is the smoothest one under given measurement data. Such solution can be obtained by performing a nonlinear constrained optimization whose Lagrangian function is of the form Z L2 ð f ; K1 Þ ¼ 0

tf

f VðtÞ2 dt þ K1 ðS −tf r2 Þ

ð9Þ

where σ is the standard deviation of the temperature reading Y(t) and Λ1 is the Lagrange multiplier that implies the reciprocal of the regularization parameter. This approach enables us to avoid the difficulty inherent in selection of the regularization of parameter. Furthermore, the first term on the right hand side implies the first order regularization and the second term is the nonlinear equality constraint that imposes the statistical consistency pffiffiffiffiffiffiffiffiffi between the measured and computed temperatures. The violation level of this constraint is given by r2 tf 2=N [1]. After accomplishing the second phase optimization, then F can be evaluated with acceptable accuracy for the third phase. The third phase starts from the result of the second phase. The second phase provides the initial guess for the optimization problem of the second phase and the prerequisite value F. The Lagrangian function of the second phase optimization, which achieves the maximum entropy solution, is expressed as Z L3 ð f ; K1 ; K2 Þ ¼ −

tf 0

 Z tf  fðtÞ 2 dt þ K1 ðS −tf r Þ þ K2 F − f ðtÞln f ðtÞdt F 0

ð10Þ

where, 0 b f (t) ≤ F, Λ1 and Λ2 are the Lagrange multiplier of each constraint. The first term on the right hand side is referred to as the information negentropy, which is opposite in sign with the information entropy [6]. The optimization in this phase aims to find the unknown temperature f (t) that maximizes the information entropy (minimizes the information negentropy) under the given temperature readings. The information entropy is analogous to the thermodynamic entropy in its mathematical properties. Generally, the entropy in an isolated system never decreases, i.e., spontaneously tends towards the maximum value possible. In other words, since MEM reflects such property of the nature, it searches the most natural solution rather than the smoothest solution. It is well known that MEM can accomplish the optimal regularization and the resolution enhancement at the same time. As is evident from the Lagrangian function, the shortcoming of this approach is that the unknown variable should be positive-definite. However, unknowns in many IHCP applications do not change sign in the entire time domain. Additionally, the optimization problem in the third phase is highly nonlinear owing to the inclusion of logarithmic functions. Thus, it demands large computational resources unavoidably. 3. Implementation The described method can be divided into two parts. The first part considers the numerical solutions of the partial differential equations involved in the formulation. The second part plays the role of the optimization utilizing the numerical solutions. Separately implementing each part and combining them completes the whole implementation. Each part is described below. All the phases requires the numerical evaluations of T(x, t; f ) and λ(x, t; f ), and the first phase calls for ΔT(x, t; f ). This work adopts finite volume method with Crank–Nicolson scheme for the discretization in time and space. The spatial domain shown in Fig. 1 is discretized with M nodal points with equal spacing. In IHCP, the time step size for the discretization usually is identical to the measurement interval Δt. Besides, all the integrals involved in this formalism are evaluated by the trapezoidal rule for accuracy. The implementation of the first phase optimization is straightforward. On the other hand, the optimization problems in the second and third phases are not easy to handle. The successive quadratic programming (SQP) is utilized for the second and third phases. SQP searches the desired optimum by repeating the following three steps consecutively: (i)

42

S.K. Kim et al. / International Communications in Heat and Mass Transfer 34 (2007) 37–44

Fig. 2. Flow chart of the three-phase MEM.

approximation of the Lagrangian function in the quadratic form and the constraints in the linear form; (ii) determination of the search direction; (iii) control of the step length along the search direction with a proper penalty. SQP requires evaluation of the Hessian matrix and the gradient vector of the objective function and the constraints. The gradient vector of the residual norm S can be obtained by solving the adjoint problem. However, the evaluation of the Hessian matrix of S is not an easy task (especially for the adjoint formulation). Sometimes, it is impractical or even impossible. Thus, in such cases, the Hessian matrix of S is approximated by the Broyden–Fletcher–Goldfarb–Shanno (BFGS) update formula [9]. Additionally, it is noted that the second constraint in the second phase is a linear constraint. In this study, SQP is implemented using an existing computer code, CFSQP [17]. This code is slightly modified to allow the user defined Hessian matrix for the objective function (CFSQP does not allow it as is). The iteration stoppage and violation level of the equality constraints are handled automatically by CFSQP. The numerical procedures are briefly explained below. As was mentioned, the entire calculation involves the three phases. In the first phase, by the CGM, the start value for the second phase is calculated. In the second phase, by the QR with the SQP, the total sum and start value for the third phase is calculated. In the third phase, the final solution is estimated using the MEM. Fig. 2 shows the flow of the tasks to be done during these three phases. The numerical discretization of the heat conduction equation is not repeated here since it is found in many literatures [1,3,15]. Furthermore, the detail about CGM is also not described here. Please, refer to Refs. [2,16]. The important thing is the total sum is calculated from the result of the second phase as shown in Fig. 2. 4. Results and discussion Test cases are investigated to verify the performance of the proposed method. The described IHCP is investigated for different forms of unknown temperature, namely rectangular and triangular forms. The unknown temperature estimated with MEM (the third phase) is compared with the original temperature and the temperatures estimated with the conjugate gradient method (the first phase) and the quadratic regularization (the second phase). Test cases have been investigated for different error levels of measurement data (σ = 0 and σ = 0.01), which is artificially generated from numerical calculations. The errors embedded to the exact data are normally distributed random data with zero mean. Besides, the tests have been performed for different time steps (Δt = 0.01 and Δt = 0.1) to examine the effect of it. The proposed method is tested for a rectangular profile including an impulse and several rising and falling edges. We have tested four different methods for the same temperature profile. They are the MEM with the QR, the MEM with the CGM, the QR and the CGM. In the MEM with the QR, the MEM estimation is conducted based on the total sum from the QR. The major contribution of this study summarized as the following two things. First, the accuracy of the total sum, F, for MEM is improved by involving the QR phase. Second, the unknown is temperature instead of heat flux [14].

S.K. Kim et al. / International Communications in Heat and Mass Transfer 34 (2007) 37–44

43

Fig. 3. Estimation of rectangular temperature profile. (a) Δt = 0.1, σ = 0, (b) Δt = 0.1, σ = 0.1, (c) Δt = 0.01, σ = 0, and (d) Δt = 0.01,σ = 0.01.

Fig. 3(a) shows the original and estimated temperature histories for Δt = 0.1 and σ = 0. The following observations can be made: (i) All the estimations are acceptable as inverse solutions; (ii) the quadratic regularization (the second phase) renders the overshoots near the edges mild; (iii) MEM enhances the resolution remarkably especially for the reconstruction of the sharp shapes such as impulse or edge. These are already found in the previous studies [13,14,18]. The MEMs with the QR and the CGM show an almost identical pattern except in the region near the final time. This is mainly owing to the inaccuracy of the total sum obtained by the CGM. The MEM with the CGM shows apparent deviation from the exact temperature near the final time. It is found that the MEM with the QR improves the accuracy of the total sum and estimation near the final time. It is already found in a previous study that the inaccuracy in the total sum is concentrated in the final time since the unknowns near the final time least affects the temperature at the sensor locations [14]. Fig. 3(b) shows results for Δt = 0.1 and σ = 0.01. We have obtained a result similar with Fig. 3(a). The following observations can be made: (i) the solutions by MEMs are smoother than those in Fig. 3(a) thanks to the regularization effect by the measurement error since making the convergence criterion smaller will induces more fluctuations and overshoots; (ii) the inaccuracy near the final time in the result by the MEM with the CGM is quite noticeable in comparison with the result by the MEM with the QR. Fig. 3(c) shows results for Δt = 0.01 and σ = 0. It reveals that: (i) CGM does not hold enough resolving power to cope with this case owing to the viscous nature of the algorithm; (ii) the quadratic regularization done by SQP, which is a kind of the Newton's method, gives more resolving power than CGM does; (iii) MEM with the QR and CGM reduces the deterministic biases near the location where the temperature changes abruptly; (iv) near the final time the MEM with the QR provides a solution more close to the exact than the MEM with the CGM. Fig. 3(d) shows results for Δt = 0.01 and σ = 0.01, which poses a quite stringent case. As can be seen in the figure, it is found that: (i) only MEM with QR and CGM can detect the first impulse; (ii) the both MEMs over-resolves the large block located in the middle owing to the so-called entropy concentration [18]; (iii) the object located in 0.295 b t b 0.355 is very difficult to recover with

44

S.K. Kim et al. / International Communications in Heat and Mass Transfer 34 (2007) 37–44

impaired temperature readings since available data is not enough to resolve it; (iv) however, MEM with the QR better resolves the object thanks to the more accurate F. In summary, the maximum entropy estimation, overall, shows good resolving power. The estimation of the total sum becomes accurate with the QR instead of the CGM. As a result, we could improve the accuracy near the final time with it.

5. Conclusion The maximum entropy method (MEM) is applied to estimation of surface temperature from temperature measurement. The inverse heat conduction problem is reformulated for MEM and a three-phase solution method utilizing the successive quadratic programming (SQP) is proposed. The inverse estimation is remarkably improved by MEM for a rectangular profile compared with the result by CGM. In this study, the total sum was calculated using both the quadratic regularization (QR) by the constrained optimization and the conventional conjugate gradient method. From the comparison of the results by theses two methods, the QR calculates the total sum more accurately. Consequently, the accuracy of the ME estimation near the final time has been improved. Acknowledgement This work was supported by Korea Research Foundation Grant (KRF-2005-003-D00039). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18]

J.V. Beck, B. Blackwell, C.R. St. Clair Jr., Inverse Heat Conduction, Wiley, 1985. A.N. Tikhonov, Y. Arsenin, Solution of Ill-posed Problems, Winstion and Sons, 1977. N. D'Souza, ASME Paper No.75-WA/HT-81, 1975. O.M. Alifanov, Inverse Heat Transfer Problem, Springer–Verlag, 1994. J.V. Beck, B. Blackwell, A. HajiSheikh, ASME Journal of Heat Transfer 39 (1996) 3649. N. Wu, The Maximum Entropy Method, Springer–Verlag, 1997. C.E. Shannon, Bell System Technical Journal 27 (1948) 379. E.T Jaynes, Physics Review 106 (1957) 620. G.V. Reklaitis, A. Ravindran, K.M. Ragsdell, Engineering Optimization, Wiley, 1983. F.M. Ramos, A. Giovannini, International Journal of Heat and Mass Transfer 38 (1995) 101. P. Lair, J. Dumoulin, P. Millan, Numerical Heat Transfer. Part A, Applications 33 (1998) 267. W.B. Muniz, F.M. Ramos, H.F.D. Velho, Computers & Mathematics with Applications 40 (2000) 1071. J. Skilling, R.K. Bryan, Monthly Notices of the Royal Astronomical Society 211 (1984) 111. S.K. Kim, W.I. Lee, International Journal of Heat and Mass Transfer 45 (2002) 381. M. Raynaud, J.V. Beck, ASME Journal of Heat Transfer 110 (1988) 30. Y. Jarny, M.N. Özişik, J.P. Bardon, International Journal of Heat and Mass Transfer 34 (1991) 2911. C. Lawrence, Z.L. Zhou, A.L. Tits, User's Guide for CFSQP Version 2.5, University of Maryland, 1997. D.L. Donoho, I.M. Johnstone, J.C. Hoch, A.S. Stern, Journal of the Royal Statistical Society. B 54 (1992) 41.