- Email: [email protected]
Estimation of thermophysical parameters revisited from the point of view of nonlinear regression with random parameters
International Journal of Heat and Mass Transfer 106 (2017) 135–141
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Estimation of thermophysical parameters revisited from the point of view of nonlinear regression with random parameters q Daniela Jarušková ⇑, Anna Kucˇerová Czech Technical University in Prague, Faculty of Civil Engineering, Thákurova 7, CZ – 166 29 Praha 6, Czech Republic
a r t i c l e
i n f o
Article history: Received 5 May 2016 Received in revised form 22 September 2016 Accepted 11 October 2016
Keywords: Estimation of thermal conductivities Non-linear regression with common factors Optimal design of experiment Monte Carlo simulations
a b s t r a c t The paper deals with the design of an experiment for estimating thermophysical parameters in a problem described by Sawaf et al. (1995) and Ruffio et al. (2012). Its aim is to illustrate potential problems connected to the application of a nonlinear regression with ‘‘nuisance” random parameters for estimating parameters of interest. Presenting four simplified versions of the problem where a solution is found numerically and its quality is assessed by a Monte Carlo simulations this paper presents typical features of nonlinear regression with random common factors. Ó 2016 Elsevier Ltd. All rights reserved.
1. Introduction
2. Experiment description
In the thermal sciences, we are often able to describe an experiment by a mathematical model. If the model presents a complete description of a physical system, it can predict the outcome of an experiment in relation to selected input parameters, e.g., thermal properties of a material. The inverse problem consists of using actual results from measurements to infer values of the input parameters [7]. As the measurements are often subject to measurement errors, we have to estimate input parameters using statistical methods. In case we are able to express a relationship between input parameters and an output by a known function f and errors of measurement are zero mean independent identically distributed (i.i.d.) random variables, we may use the method of least squares to estimate the input parameters. The output often depends not only on the input parameters, but also on some additional parameter(s) that may be chosen before an experiment starts. We call them design parameters. The goal of stochastic analysis is to choose the design parameters for estimating the input parameters as precisely as possible. The paper examines design experiment problems in nonlinear regression.
Our paper was inspired by Ruffio et al. [4] who presents an analysis of a numerical experiment designed by Sawaf et al. [5]. The aim of the experiment was to identify the thermal parameters of an orthotropic homogeneous material - more precisely, a constant thermal conductivity along the x-axis kx , a thermal conductivity along the y-axis ky , and a thermal volumetric capacity C. The experiment can be described as follows. A square sample is exposed to a constant and uniform heat flux u on the left and bottom boundaries, while the right and top edges are insulated. The system obeys the following equation:
q This work was supported by the Czech Science Foundation under grants 15-09663S and 15-07299S. ⇑ Corresponding author. E-mail addresses: [email protected] (D. Jarušková), Anna.Kucerova@fsv. cvut.cz (A. Kucˇerová).
http://dx.doi.org/10.1016/j.ijheatmasstransfer.2016.10.044 0017-9310/Ó 2016 Elsevier Ltd. All rights reserved.
C
@T @2T @2T ¼ kx 2 þ ky 2 ; @t @x @y
0 6 x 6 Lx ; 0 6 y 6 Ly ; 0 6 t 6 s:
ð1Þ
The boundary and initial conditions are defined by:
kx @T ðx ¼ 0Þ ¼ u; @x
ky @T ðy ¼ 0Þ ¼ u; @y
ðx ¼ Lx Þ ¼ 0; ky @T ðy ¼ Ly Þ ¼ 0; kx @T @x @y Tðx; y; 0Þ ¼ 0: It is assumed that temperature is measured by one or more sensors at equidistant time points. Ruffio et al. [4] as well as Sawaf et al. [5] consider an experiment with measurements taken at times t ¼ 1; . . . ; 60 (s). The number of sensors was fixed (in both papers, three sensors are used). The aim of an analysis of the numerical model (1) is to determine position(s) of sensor(s) to identify thermal properties of the material as accurately as possible.
D. Jarušková, A. Kucˇerová / International Journal of Heat and Mass Transfer 106 (2017) 135–141
136
Nomenclature thermal conductivity along the xaxis thermal conductivity along the yaxis true value of kx true value of ky thermal volumetric capacity temperature field sample dimensions u heat flux density ru standard deviation of u J number of sensors x x-coordinate of a sensor position y y-coordinate of a sensor position DxðjÞ; j ¼ 1; . . . ; J random displacements along the x-axis DyðjÞ; j ¼ 1; . . . ; J random displacements along the y-axis rx standard deviation of DxðjÞ ry standard deviation of DyðjÞ n number of measurements t i ; i ¼ 1; . . . ; n time points of measurement
kx ky kx ky C T Lx ; Ly
Ruffio et al. [4] claim that a solution of (1) may be found analytically in the form of a quickly converging series. In our paper, we assume that a final solution may be quite accurately approximated by the first term only, i.e.,
Tðx; y; tÞ ¼ hx ðt; kx ; C; u; xÞ þ hy ðt; ky ; C; u; yÞ;
where
u pffiffi p hx ðt; kx ; C; u; xÞ ¼ p2ffiffiffiffiffi t F exffit ; Ckx e u ffi pffiffi p t F yffi hy ðt; ky ; C; u; yÞ ¼ p2ffiffiffiffiffi t
Cky
pffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffi with e x ¼ ðx=2Þ C=kx , e y ¼ ðy=2Þ C=ky and
FðzÞ ¼
Z expðz2 Þ 2 pffiffiffiffi z 1 pffiffiffiffi
p
p
0
z
2 ev dv ;
z P 0:
Because sensors measure temperature with measurement errors, we may assume that the behavior of a measured temperature can be modelled using a nonlinear regression and the problem described above is a design experiment problem in the nonlinear regression [6]. In the case of three sensors, it was numerically solved in [5]. However, Ruffio et al. [4] argue that measured temperature is not only affected by measurement errors but other input parameters may be random as well. (The very same problem was also considered by Emery et al. [1]). For instance, heat flux density applied to sample boundaries is not perfectly known. Similarly, a sensor may not be located exactly at a designed position. Then, instead of being measured at a point ðx; yÞ temperature is measured at a point ðx þ Dx; y þ DyÞ, where Dx and Dy are random. Ruffio et al. [4] assume that the distribution of the variables u; DxðjÞ; DyðjÞ; j ¼ 1; . . . ; J (J being a number of sensors), is known. In this way, they deal with the design of an experiment in a nonlinear regression with random parameters. Additionally, it is implicitly assumed that random parameters do not change during one experiment. In the other words, each of them attains only one fixed value for all time points t ¼ 1; . . . ; 60, but they may attain different values if the experiment is repeated. This means that these random parameters play the role of so-called random common factors. Then, as we illustrate below, the variability of the estimates of the parameters of interest is practically entirely determined by the variability of random parameters, i.e. it depends on the variability of measurement errors only very slightly. The goal of an analysis is to find the ‘‘optimal” position(s) of sensor(s) so that if
s final time of measurement Y i ; i ¼ 1; . . . ; n values of the dependent variable b vector of the parameters of interest b true value of b c vector of random parameters c0 mean of c Rc variance–covariance matrix of c n design parameter rm standard deviation of measurement errors f ; g; h regression functions F matrix of partial derivatives of f with respect to the components of b computed at values b and c0 Z matrix of partial derivatives of f with respect to the components of c computed at values b and c0 v21a ðpÞ ð1 aÞ 100% quantile of the v2 distribution with p degrees of freedom MSEðbÞ mean square error
we repeat the experiment, the differences between estimates of parameters and their true values are small. If we planned to do only one experiment, it would be better to consider the random parameters to be nuisance parameters that might be estimated. Their known distribution might serve as a prior distribution and we might try to find the estimates of the parameters of interest as well as of the nuisance parameters using a Bayesian approach. This is, however, a different statistical problem than we are solving here. Despite general agreement with the method for designing an experiment as suggested by Ruffio et al. [4], we are afraid that the authors do not fully realize problems connected to nonlinear regression models with random factors. Our objections do not only concern the specific problem described above, but all similar problems. For the sake of clarity, we simplified the situation assuming that the parameters of interest are heat conductivities kx and ky only, while thermal volume capacity C is known to be equal to 1,700,000 (Jm3 K1 ). 3. Non-linear regression with non-random parameters We suppose that we observe a sequence of pairs fðt i ; Y i Þ; i ¼ 1; . . . ; ng. We assume that ft i ¼ is=n; i ¼ 1; . . . ; ng are time points from the interval ½0; s and
Y i ¼ f ðt i ; bÞ þ ei ;
i ¼ 1; . . . ; n;
ð2Þ
where f ðt; bÞ is a known regression function of t 2 ½0; s and a p-dimensional vector of unknown parameters of interest b ¼ ðb1 ; . . . ; bp ÞT . The random errors fei g are i.i.d. random variables with a variance r2m . We suppose that the true value of b is b 2 ½bl ; bu and we estimate b using the least squares method: n X 2 b b ¼ argmin ðY i f ðt i ; bÞÞ : b2½bl ;bu i¼1
Under a set of general conditions, b b is a consistent estimate of b pffiffiffi P (i.e. b b ! b as n ! 1) and nð b b b Þ has asymptotically a 1 p dimensional normal distribution N 0; ðF T F Þ r2m , see [6]. The matrix F is a matrix of partial derivatives of f with respect to b at the right value b .
D. Jarušková, A. Kucˇerová / International Journal of Heat and Mass Transfer 106 (2017) 135–141
Further, we assume that the regression function f ðt; b; nÞ depends on some parameter n. The goal of our analysis is to choose n in such a way that the convergence of b b to b is as fast as possible. As the volume of a ð1 aÞ100% asymptotic confidence ellipsoid
T ðb b b Þ F T F ð b b b Þ 6 r2m v21a ðpÞ is proportional to
of b b under an assumption that at the same time the parameter c has been randomly chosen from a distribution N c0 ; r2c and the measurement errors fei ; i ¼ 1; . . . ; ng have been randomly chosen b by from Nð0; r2 Þ. Here, we can express b m
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 det F T F . The D criterion consists in mini-
mizing a determinant of an asymptotic variance–covariance matrix 1 . We will of b b and that is here the same as to minimize det F T F also consider the F criterion, introduced by Ruffio et al. [4], where P 1=2 p 2 b is to be minimized. i¼1 Var b i =ðbi Þ 1 The asymptotic variance–covariance matrix F T F r2m depends on the right value b so that an optimal value nop with respect to the D as well as the F criterion depends on b as well. This is a substantially different than in linear regression case where a choice of nop is independent of b , see [3]. One way to overcome the problem is to proceed in two steps, i.e. in the first step, estimating the parameter b only roughly and in the second step, finding nop with respect to the estimate obtained in the first step. If we know that b 2 ½bl ; bu , we can minimize a criterion with respect to the ‘‘worse” b in the interval ½bl ; bu . This is basically the idea of Ruffio et al. [4]. 4. Non-linear regression with random parameters
P is is P is h n ;n g n ;n g ;n e b P 2 is b ¼ b þ c c0 þ P 2nis i : g n ;n g n ;n
As c and fei g are independent, the unconditional distribution of b b is normal with the mean b and a variance:
P is is !2 h n ;n g n ;n r2 P 2 is r2c þ P 2 mis : g n ;n g n ;n
Var b b¼
i ¼ 1; . . . ; n: ð3Þ
: b¼ Var b
Rs 0
hðt; nÞgðt; nÞdt Rs g 2 ðt; nÞdt 0
!2
r2c þ
r2m n
Rs 0
1 : g 2 ðt; nÞdt
may minimize either (4) or (5). For a large n, the difference will be small. In the case of a nonlinear regression model with a random parameter c:
i ¼ 1; . . . ; n;
ð6Þ
Ruffio et al. [4] suggest replacing the model (6) with a linear model
Y i ¼ f ðt i ; b ; c0 ; nÞ þ þ
@f ðt i ; b ; c0 ; nÞðb b Þ @b
@f ðt i ; b ; c0 ; nÞðc c0 Þ þ ei ; @c
i ¼ 1; . . . ; n
@f ðt ; b ; @c i
parameter c is independent of random errors fei ; i ¼ 1; . . . ; ng that are i.i.d. random variables. For simplicity, we suppose that they are distributed according to Nð0; r2m Þ. Suppose that in our experiment, the parameter c actually c that might be different from c0 . As we do not attained a value e Pn know this value, we estimate b by minimizing i¼1 ðY i l
1 1 1 F T F F T Z Rc Z T F F T F þ r2m F T F ;
b b ¼ b þ
Pn
is i¼1 ðY i Þg n ; n Pn 2 is ;n i¼1 g n
l
:
The conditional distribution of b b given c ¼ e c describes a stochastic behavior of b b under the assumption that the experiment is perc . It is normal with formed many times under a condition that c ¼ e a mean
b þ
P is is h n ;n g n ;n e P 2 is c c0 g n ;n
@f This is indeed model (3) with gðt i ; nÞ ¼ @b ðt i ; b ; c0 ; nÞ and hðt i ; nÞ ¼
c0 ; nÞ.
For a multidimensional parameter b, we obtain an approximate b in the form presented by Ruffio variance–covariance matrix of b et al. [4]:
r2m
g 2 ins ; n
:
The unconditional distribution of b b under an assumption that c has a normal distribution N c0 ; r2c describes a stochastic behavior
ð7Þ
where F is a matrix of partial derivatives of f with respect to components of b, Z is a matrix of partial derivatives of f with respect to components of c, both computed at b and c0 , and Rc is a diagonal matrix with a diagonal ðVarc1 ; . . . ; Varck Þ. Ruffio et al. [4] suggest applying either the D criterion or the F criterion to the matrix (7). Similarly as in the linear regression with random parameters, b is a non-consistent estimator. Moreover, b b is no the estimate b longer an unbiased estimator and its bias does not disappear as the number of observations tends to infinity. The reason is the b of the least squares following. We are looking for a minimum b function MSEðbÞ:
MSEðbÞ ¼
and a variance
P
ð5Þ
The estimate b b is an unbiased estimator but it is not consistent Rs unless 0 hðt; nÞgðt; nÞdt ¼ 0. The asymptotic distribution of b b b is a zero mean normal distribution with a variance equal to the first term in (5). In case we wish to minimize Var b b with respect to n, we
We suppose that gðt; nÞ and hðt; nÞ are continuous functions on ½0; s ½nl ; nu . We suppose that l is a known real parameter, while b is a parameter that is to be estimated. Further, we suppose that c is a random parameter distributed according to Nðc0 ; r2c Þ. The
gðis=n; nÞðb b ÞÞ2 . Therefore, the sought estimate has a form:
ð4Þ
If n is large, the variance Var b b behaves approximately as
Y i ¼ f ðt i ; b; c; nÞ þ ei ; In this section we explain properties of a least squares estimate of a parameter of interest in a non-linear regression with random parameters. We start analyzing a linear regression model in the form:
Y i ¼ l þ gðis=n; nÞðb b Þ þ hðis=n; nÞðc c0 Þ þ ei ;
137
¼
n 2 1X Y i f ðt i ; b; c0 ; nÞ n i¼1 n 2 1X f ðti ; b ; c; nÞ f ðt i ; b; c0 ; nÞ n i¼1
þ
n n 2X 1X ei ðf ðt i ; b ; c; nÞÞ f ðt i ; b; c0 ; nÞÞ þ e2 : n i¼1 n i¼1 i
ð8Þ
D. Jarušková, A. Kucˇerová / International Journal of Heat and Mass Transfer 106 (2017) 135–141
138
As n tends to infinity, the second term of (8) converges to zero and Rs b minimizes ðf ðt; b ; c; nÞ the third term to r2 . The optimal b m
0
f ðt; b; c ; nÞÞ dt and it is generally a nonlinear function of c. Due to b may differ from b , nonlinearity, the mean of the distribution of b 0
2
b may differ from the matrix and the variance–covariance matrix of b (7). Clearly, this discrepancy does not disappear as n tends to infinb canity. Therefore, an approximate variance–covariance matrix of b b unless not be a criterion for determining accuracy of the estimate b it is approximately a linear function of c. For a one dimensional parameter, we should compare the quality of an estimate b b using 2 a square mean deviation from the true value E ð b b b Þ . We recall 2 2 b þ ðE b that Eð b b b Þ ¼ Var b b b Þ . For a p-dimensional paramePp 2 2 b ter, we consider the F criterion, where is to j¼1 Eð bj bj Þ =bj
be minimized. However, we cannot evaluate it unless we know the true value b and we get into a vicious circle.
we infer the matrix F comprising of two columns, namely x ¼ ðx1 ; . . . ; xn ÞT and y ¼ ðy1 ; . . . ; yn ÞT , where for i ¼ 1; . . . ; n
pffiffiffiffi e 2u x xi ¼ pffiffiffi 3=2 ti G pffiffiffiffi ; ti C ðkx Þ with e x ¼ ðx=2Þ
pffiffiffiffi e 2u y yi ¼ pffiffiffi 3=2 t i G pffiffiffiffi ti C ðky Þ
qffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffi y ¼ ðy=2Þ C=ky and kx , ky being true values of C=kx , e
parameters kx , resp. ky . First, we consider the D criterion. Without a loss of generality we confine ourselves to a special case of kx ¼ ky ¼ 1 and search b y op ð1; 1Þ as for general values of the conductivities kx x op ð1; 1Þ; b and ky . It holds that
b x op ðkx ; ky Þ ¼
qffiffiffiffiffi kx b x op ð1; 1Þ;
b y op ðkx ; ky Þ ¼
qffiffiffiffiffi ky b y op ð1; 1Þ
ð10Þ
and for
1
T ; Dop ð1; 1Þ ¼ min det ðF ð1; 1ÞÞ ðF ð1; 1ÞÞ x;y
5. Numerical results
(
In this section we highlight some interesting properties of the design of an experiment in a nonlinear regression. We study model (1) with Lx ¼ Ly ¼ 0:05 (m), see Fig. 1. Recall that the unknown parameter of interest is the vector of thermal conductivities ðkx ; ky ÞT . We start with a non-linear regression where no random factor is present. Case 1 We are to find a position ðx; yÞ for one sensor to estimate b ¼ ðkx ; ky ÞT . We assume that all other parameters are fixed; the
Dop
¼ min det
x;y
F
ðkx ; ky Þ
T
F
ðkx ; ky Þ
) 1 ;
it holds
3 Dop ¼ kx ky Dop ð1; 1Þ:
ð11Þ
We proved analytically that for kx ¼ ky ¼ 1 there exist two optimal symmetric positions that are situated on the boundaries of the sample. Looking for an optimal value of the x-coordinate of the first optimal position on a grid ½0 : 0:0001 : 0:05 we have found that ð1Þ
ð1Þ
volumetric thermal capacity C ¼ 1; 700; 000 ðJm3 K1 Þ, the heat
ðxop ¼ 0:0023; yop ¼ 0Þ. Clearly, for the second optimal position it
flux density u ¼ 25; 000 ðWm2 Þ, and the sensor is located exactly at ðx; yÞ so that Dx ¼ 0 and Dy ¼ 0. In agreement with [4], temperature is measured at times t i ¼ i; i ¼ 1; . . . ; n ¼ 60. A standard deviation of measurements errors rm ¼ 0:1. Under these assumptions, the regression function’s form is:
holds
e e 2u pffiffi 2u pffiffi x y f ðt; kx ; ky ; x; yÞ ¼ pffiffiffiffiffiffiffiffi tF pffiffi þ pffiffiffiffiffiffiffiffi t F pffiffi : Ckx Cky t t
ð9Þ
For finding an optimal ðx; yÞ, we compute the asymptotic variance– covariance matrix of b b. Introducing a new function GðzÞ by
Z 2 GðzÞ ¼ FðzÞ þ z 1 pffiffiffiffi
p
0
z
2 ev dv ;
ð2Þ ðxop
ð2Þ
¼ 0; yop ¼ 0:0023Þ. The minimal value of the D criterion 13
is 1:40 10 . For kx ¼ 0:6 and ky ¼ 4:7, considered in [4], the optimal positions of a sensor and the optimal value of the D criterion are pffiffiffiffiffiffiffi ð1Þ ð1Þ ð2Þ ð2Þ indeed xop ¼ 0:6 0:0023 ¼ 0:0018; yop ¼ 0 and xop ¼ 0; yop ¼ pffiffiffiffiffiffiffi 4:7 0:0023 ¼ 0:0050, Dop ¼ ð0:6 4:7Þ3 1:40 1013 ¼ 3:14 1012 . Notice, see (10), that an optimal position may be substantially affected by the true values of parameters of interest. The criterion F is not homogeneous with respect to kx and ky and that is why we have to find the optimal positions for kx ¼ ky ¼ 1 and kx ¼ 0:6; ky ¼ 4:7 separately. For kx ¼ ky ¼ 1 the optimal posi ð1Þ ð1Þ ð2Þ ð2Þ tions are xop ¼ 0:0025; yop ¼ 0 and xop ¼ 0; yop ¼ 0:0025 with an optimal value 1:8 103 . We see that in this case and kx ¼ ky ¼ 1 the optimal positions for the criteria D and F are not the same but they differ only slightly. For kx ¼ 0:6 and ky ¼ 4:7 only one optimal solution xop ¼ 0:0024; yop ¼ 0 exists and the optimal value of the criterion F is 0:0016. For kx ¼ ky ¼ 1 the matrix ðF ð1; 1ÞÞ ðF ð1; 1ÞÞ is singular when a sensor is placed on the diagonal of the square sample, i.e. x ¼ y, as the both coordinates of F ð1; 1Þ are the same. Positions near the diagonal yield an ill-conditioned matrices, e.g. the ratio of eigenT
ð1Þ
ð1Þ
values for ðxop ¼ 0:0023; yop ¼ 0Þ is 75.8. For a general kx ; ky the qffiffiffiffiffi pffiffiffiffiffi matrix F T F is singular when x= kx ¼ y= ky . To avoid the problem with singularities instead of minimizing the D and F criterion we maximized their inverses. In the three following cases, we estimate the parameter
Fig. 1. Description of experiment to optimize.
b ¼ ðkx ; ky ÞT in the nonlinear regression model (6) with the regression function (9) and a standard deviation of measurement errors rm ¼ 0:1, but we suppose that some of the input parameters are random. We compare the value of the F criterion for the ‘‘optimal” sensor position where the variances were obtained using the lin-
D. Jarušková, A. Kucˇerová / International Journal of Heat and Mass Transfer 106 (2017) 135–141
earization with the value of the F s , resp. F s , criterion where means and variances were obtained by numerical solutions of a non-linear minimization problem with parameters obtained by the Monte Carlo simulations. Our aim is to show that in some cases ‘‘good design” obtained by linearization yields very bad design as the mean square errors of the least squares estimates of the identified parameters are large. The optimal position(s) of sensor(s) were found numerically on a grid ½0 : 0:0001 : 0:05 ½0 : 0:0001 : 0:05 or using an evolutionary GRADE algorithm extended by niching strategy CERAF [2]. The matrices used in the D and F criteria were the approximate variance–covariance matrices given in (7). After finding the optimal sensor position, we simulated 100,000 nonlinear models (6) with random common factors. Using the least squares method, we estimated the parameter of interest b for all 100,000 simulated models. For obtaining the least squares estib we applied the Matlab function ‘‘lsqcurvefit” with minimates b, mum precision of the target least squares function ‘TolFun’ equal to 1025 and minimum precision of the estimated parameters ‘TolX’ equal to 1015 . Case 2 We consider the nonlinear regression model (6) with the regression function (9) and the true values of parameters kx ¼ ky ¼ 1. We suppose that heat flux density u is a random variable distributed according to a normal distribution with a mean
u0 ¼ 25; 000 and ru ¼ 50 ðWm2 Þ. All other parameters are fixed.
We are looking for the position of one sensor ðx; yÞ. First, we consider the D criterion. Fig. 2 shows contours of detðVar b bÞ on ½0; 0:005 ½0; 0:005 computed by (7) with respect
to ðx; yÞ. The optimal value is attained at two positions x ¼ 0:0027; y ¼ 0 and x ¼ 0; y ¼ 0:0027 and the attained minimum is 3:04 1011 . Using the results of simulations for x ¼ 0:0027; y ¼ 0 k y 1Þ ¼ 0:0006 and we estimate the bias Eðb k x 1Þ ¼ 0:0002; Eðb k y ¼ 0:0008. the variances Var b k x ¼ 0:0002; Var b For the F criterion, the optimal value is attained at two positions x ¼ 0:0023; y ¼ 0 and x ¼ 0; y ¼ 0:0023. The attained minimum is 0:0036. Using the results of simulations for x ¼ 0:0023; y ¼ 0, we estimate the bias Eðb k x 1Þ ¼ 0:0002; Eðb k y 1Þ ¼ 0:0006, and the k y ¼ 0:0008 and the correlation coefvariances Var b k x ¼ 0:0002; Var b b k y Þ ¼ 0:96. ficient Corrð k x ; b We have seen that the accuracy of the parameters’ estimates in regression with random parameters is given by the variability of the random factors. Because the variability of the heat flux density is small, the optimal sensor positions for a random and nonrandom u are very close. Moreover, the dependence of b k x and kby on u is practically linear and therefore the biases of the estimates are negligible. This is illustrated by Fig. 3 that shows the dependence of b k x k and kby k on u 2 ½u0 3ru ; u0 þ 3ru for x
y
x ¼ 0:0023; y ¼ 0; kx ¼ 1 and ky ¼ 1. −3
5
y.10
Fig. 3. Dependence of the estimates b k x kx (solid line) and b k y ky (dashed line) on ufor kx ¼ ky ¼ 1 and x ¼ 0:0023; y ¼ 0.
Case 3 We consider again the nonlinear regression model (6) with the regression function (9) with the true values of parameters kx ¼ 0:6 and ky ¼ 4:7. We would like to measure temperature using one sensor and we know that we are not able to place it exactly at the designed position. The random displacements Dx and Dy from the right position ðx; yÞ are i.i.d., distributed according to a normal distribution with zero mean and a standard deviation rx ¼ ry ¼ 0:0005 (m). All other parameters are fixed. Using the F criterion, the optimal solution is attained at x ¼ 0:0165; y ¼ 0:0158 and the minimum of the criterion F is 0.0387. The contours of the F criterion are presented in Fig. 4. Notice that the ‘‘optimal” position differs substantially from the optimal position when the random common factors are not present. Using simulated values, the biases estimates are Eðb k x 0:6Þ ¼ 0:0001 and Eðb k y 4:7Þ ¼ 0:0436. The variances are k y ¼ 0:2486 and Corrðb kx ; b k y Þ ¼ 0:87. The value Var b k x ¼ 0:0117; Var b of the F s criterion is equal to 0.2093 and that of F s to 0.2095. We see that the biases are small but the values of F s , resp. of F s criterion, obtained by simulations, differ from the values of the F criterion obtained using (7). Finally, we simulated model (6) without measurement errors, i.e. r2m ¼ 0. For x ¼ 0:0165 and y ¼ 0:0158, the estimates of the biases and the variances changed only very little. The biases are 0.0000 and 0.0403; the variances are 0.0111 and 0.2483. We see that if n is large, the random measurement errors do not affect variability of the parameters’ estimates. Case 4 In this case we compare the quality of four designs suggested by Ruffio et al. [4] and denoted by U1; U2; U3, and U4. All four designs correspond to the location of three sensors. The aim of the analysis presented in [4] was to find the optimal positions not only to estimate thermal conductivities but also the volumetric thermal capacity C. In accordance with the previous cases, we assess the designs according to a quality of estimates of thermal conductivities kx and ky assuming C ¼ 1; 700; 000 Jm3 K1 . Further, in agreement with [4] we assume that u is a normally distributed random y [m] 0.03
[m]
4
0.025
3
0.02
2
0.015
1 0
139
0
1
2
3
4 5 −3 x.10 [m]
Fig. 2. Contours of det Var b b computed using an approximate variance–covariance matrix (7) for ðx; yÞ 2 ½0; 5 103 ½0; 5 103 .
0.01 0.01
0.015
0.02
0.025
0.03 x [m]
Fig. 4. Contours of criterion F computed using an approximate variance–covariance matrix (7) for ðx; yÞ 2 ½0:01; 0:03 ½0:01; 0:03.
D. Jarušková, A. Kucˇerová / International Journal of Heat and Mass Transfer 106 (2017) 135–141
140
variable with a mean u ¼ 25; 000 and a standard deviation
ru ¼ 50 (W m2 Þ, and that, for j ¼ 1; . . . ; 3, the variables DxðjÞ as
well as DyðjÞ had zero mean normal distribution with a standard deviation of 0.0005. All random variables are independent. The true value of kx is 0.6 and of ky is 4.7. Paper [4] does not provide the exact positions of the sensors for U1; U2; U3, and U4 and we interfered them from Fig. 10, [4]. We assume the positions to be: U1
U3
x1 ¼ 0:011 x2 ¼ 0:011 x3 ¼ 0:037
y1 ¼ 0:031 y2 ¼ 0:031 y3 ¼ 0:05
U2
x1 ¼ 0:010 x2 ¼ 0:033 x3 ¼ 0:044
y1 ¼ 0:031 y2 ¼ 0 y3 ¼ 0
U4
x1 ¼ 0:033 x2 ¼ 0:044 x3 ¼ 0:016
y1 ¼ 0:023 y2 ¼ 0:023 y3 ¼ 0:05
x1 ¼ 0:011 x2 ¼ 0:030 x3 ¼ 0:039
y1 ¼ 0:048 y2 ¼ 0 y3 ¼ 0
We can see that all designs have either one or two sensors on the lower or upper boundary of a sample. In such a case, the corresponding Dy cannot have a zero mean normal distribution. Therefore, we assume that if a sensor was located on the lower boundary y ¼ 0, then Dy had a (positive) half normal distribution with a standard deviation stdðDyÞ ¼ 0:0005. If it was located on the upper boundary, then Dy had a (negative) half normal distribution with a standard deviation stdðDyÞ ¼ 0:0005. In such a way, the mean of ycoordinate of a sensor originally located on the lower boundary is 6:6 104 , and of a sensor originally located on the upper boundary is 0:05 6:6 104 . Table 1 presents values of F based on the approximate variances obtained by (7) and the values of F s and F s obtained by simulations. All biases are small. However, by ordering the designs according to the performance of F based on linearization, we get: U2; U1; U4, and U3 while by ordering according to the performance of F s , resp. F s , we get: U2; U4; U1; U3. This is caused by the difference between the values of F and F s . Using a genetic GRADE algorithm extended by niching strategy CERAF [2] we found an ”optimal” design U5 with respect to the F criterion:
x1 ¼ 0:025327;
y1 ¼ 0:024367;
x2 ¼ 0:002705;
y2 ¼ 0:012212;
x3 ¼ 0:004806;
y3 ¼ 0:005897:
Using (7) the value of the F criterion is equal to 0.06 as Var b k x ¼ 0:0006 and Var b k y ¼ 0:0449. Using simulations we estib k y equal to mated the bias of k x equal to 0.129 and the bias of b k y ¼ 13:13; F s ¼ 0:973 1.078. Moreover, we get Var b k x ¼ 0:127; Var b
and F s ¼ 1:019. Clearly, the design U5 is a bad design comparing it with U1; U2; U3, and U4 despite the value of the F criterion computed by (7) being very small. Not only the estimates have extremely large biases but the variances are large as well. Fig. 5 illustrates nonlinearity of relationship between b k x , resp. b k y , and
the x-coordinate of the third sensor. We would like to remark that for all analyzed designs a difference between the corresponding elements of the variance–
Table 1 The values of F based on the approximate variances obtained by (7) and the values of F s and F s obtained by simulations. U1 U2 U3 U4
F F F F
¼ 0:1285 ¼ 0:0964 ¼ 0:1649 ¼ 0:1478
Fs Fs Fs Fs
¼ 0:1470 ¼ 0:0877 ¼ 0:1549 ¼ 0:1400
F s F s F s F s
¼ 0:1473 ¼ 0:0883 ¼ 0:1546 ¼ 0:1404
Fig. 5. Dependence of the estimates ðb k x kx Þ=kx (solid line) and ðb k y ky Þ=ky (dashed line) on x3 for kx ¼ 0:6; ky ¼ 0:8 and the design U5.
covariance matrix for the model with and without measurement errors is less 103 . Finally, we study properties of the robust version of the criterion F. If it is known that the parameters of interest are elements of some known intervals, i.e.,
b 2 ½b1 db1 ; b1 þ db1 . . . ½bp dbp ; bp þ dbp ; the authors [4] suggest evaluating the F criterion not only for b ¼ ðb1 ; . . . ; bp ÞT but also for fðb1 db1 ; . . . ; bp dbp ÞT g and compute the maximum of all values of F as a worst case scenario. We supposed that kx 2 ½0:3; 0:7 and ky 2 ½3:0; 7:0 and computed the values of F s and F s obtained by simulations for all designs U1; U2; U3, and U4 and ðkx ; ky Þ ¼ ð0:6; 4:7Þ; ð0:3; 3:0Þ; ð0:3; 7:0Þ; ð0:7; 3:0Þ; ð0:7; 7:0Þ. The designs U1; U3 and U4 are relatively unsensitive to the change of true values of kx and ky , but the design U2 is sensitive to the change of kx . Applying a genetic GRADE algorithm [2] we have found the ‘‘optimal design” U6 with the respect to FR criterion:
x1 ¼ 0:040665;
y1 ¼ 0:034156;
x2 ¼ 0:014721;
y2 ¼ 0:046086;
x3 ¼ 0:014148;
y3 ¼ 0:05:
The value of F is 0.089, and of FR is 0.091. Using simulations we found that the biases for kx 2 ½0:3; 0:7 and ky 2 ½3:0; 7:0 are generally small. The value of F s is 0.073 and of its robust version is 0.160. Therefore, we consider the U6 design to be very good. 6. Conclusion We numerically analyzed a heat transfer problem described by Sawaf et al. [5] and Ruffio et al. [4] with the aim of finding optimal sensor positions for estimating thermal conductivities. We have shown that if we use the approach suggested in [4] with a criterion based on an approximate variance–covariance matrix of estimates in nonlinear regression models with random parameters, we encounter some problems. First, if the number of measurements n is large, the level of inaccuracy of estimated parameters is almost completely determined by the variability of random parameters (i.e. common factors). Moreover, the estimates of the parameters of interest are not consistent, i.e. the estimate does not tend to the true value as the number of measurements increases. In such a case, their possible bias that arises from a nonlinear dependence of the estimates of parameters of interest on random ‘‘nuisance” parameters may play an important role. Non-linearity may also cause the substantial difference between the variance–covariance matrix obtained using linearization from the true variance–covariance matrix of the estimates. In some problems, approximate linearity may be analyzed by repeating numerical calculations for those parameters
D. Jarušková, A. Kucˇerová / International Journal of Heat and Mass Transfer 106 (2017) 135–141
values that are attained most ‘‘likely”. An alternative method is to assess the quality of estimates using Monte Carlo simulations. It would be difficult to find an optimal design with simulations, but at least the characteristics of the solution obtained using linearization should be compared with the characteristics obtained via simulations. References [1] A.F. Emery, A.V. Nenarokomov, T.D. Fadale, Uncertainties in parameter estimation: the optimal experiment design, Int. Heat Mass Transfer 43 (18) (2000) 3331–3339.
141
[2] A. Ibrahimbegovic, C. Knopf-Lenoir, A. Kucˇerová, P. Villon, Optimal design and optimal control of structures undergoing finite rotations and elastic deformations, Int. J. Numer. Meth. Eng. 61 (2004) 2428–2460. [3] L. Pronzato, A. Pázman, Design of Experiments in Nonlinear Models, Springer, New York, 2013. [4] E. Ruffio, D. Saury, D. Petit, Robust experiment design for the estimation of thermophysical parameters using stochastic algorithms, Int. Heat Mass Transfer 55 (2012) 2901–2915. [5] B. Sawaf, M.N. Ozisik, Y. Jarny, An inverse analysis to estimate linearly temperature dependent thermal conductivity components and heat capacity of an orthotropic medium, Int. Heat Mass Transfer 38 (16) (1995) 3005–3010. [6] G.A.F. Seber, C.J. Wild, Nonlinear Regression, Wiley, New York, 1989. [7] A. Tarantola, Inverse Problem Theory and Methods for Model Parameter Estimation, Siam, Philadelphia, 2005.
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Estimation of thermophysical parameters revisited from the point of view of nonlinear regression with random parameters q Daniela Jarušková ⇑, Anna Kucˇerová Czech Technical University in Prague, Faculty of Civil Engineering, Thákurova 7, CZ – 166 29 Praha 6, Czech Republic
a r t i c l e
i n f o
Article history: Received 5 May 2016 Received in revised form 22 September 2016 Accepted 11 October 2016
Keywords: Estimation of thermal conductivities Non-linear regression with common factors Optimal design of experiment Monte Carlo simulations
a b s t r a c t The paper deals with the design of an experiment for estimating thermophysical parameters in a problem described by Sawaf et al. (1995) and Ruffio et al. (2012). Its aim is to illustrate potential problems connected to the application of a nonlinear regression with ‘‘nuisance” random parameters for estimating parameters of interest. Presenting four simplified versions of the problem where a solution is found numerically and its quality is assessed by a Monte Carlo simulations this paper presents typical features of nonlinear regression with random common factors. Ó 2016 Elsevier Ltd. All rights reserved.
1. Introduction
2. Experiment description
In the thermal sciences, we are often able to describe an experiment by a mathematical model. If the model presents a complete description of a physical system, it can predict the outcome of an experiment in relation to selected input parameters, e.g., thermal properties of a material. The inverse problem consists of using actual results from measurements to infer values of the input parameters [7]. As the measurements are often subject to measurement errors, we have to estimate input parameters using statistical methods. In case we are able to express a relationship between input parameters and an output by a known function f and errors of measurement are zero mean independent identically distributed (i.i.d.) random variables, we may use the method of least squares to estimate the input parameters. The output often depends not only on the input parameters, but also on some additional parameter(s) that may be chosen before an experiment starts. We call them design parameters. The goal of stochastic analysis is to choose the design parameters for estimating the input parameters as precisely as possible. The paper examines design experiment problems in nonlinear regression.
Our paper was inspired by Ruffio et al. [4] who presents an analysis of a numerical experiment designed by Sawaf et al. [5]. The aim of the experiment was to identify the thermal parameters of an orthotropic homogeneous material - more precisely, a constant thermal conductivity along the x-axis kx , a thermal conductivity along the y-axis ky , and a thermal volumetric capacity C. The experiment can be described as follows. A square sample is exposed to a constant and uniform heat flux u on the left and bottom boundaries, while the right and top edges are insulated. The system obeys the following equation:
q This work was supported by the Czech Science Foundation under grants 15-09663S and 15-07299S. ⇑ Corresponding author. E-mail addresses: [email protected] (D. Jarušková), Anna.Kucerova@fsv. cvut.cz (A. Kucˇerová).
http://dx.doi.org/10.1016/j.ijheatmasstransfer.2016.10.044 0017-9310/Ó 2016 Elsevier Ltd. All rights reserved.
C
@T @2T @2T ¼ kx 2 þ ky 2 ; @t @x @y
0 6 x 6 Lx ; 0 6 y 6 Ly ; 0 6 t 6 s:
ð1Þ
The boundary and initial conditions are defined by:
kx @T ðx ¼ 0Þ ¼ u; @x
ky @T ðy ¼ 0Þ ¼ u; @y
ðx ¼ Lx Þ ¼ 0; ky @T ðy ¼ Ly Þ ¼ 0; kx @T @x @y Tðx; y; 0Þ ¼ 0: It is assumed that temperature is measured by one or more sensors at equidistant time points. Ruffio et al. [4] as well as Sawaf et al. [5] consider an experiment with measurements taken at times t ¼ 1; . . . ; 60 (s). The number of sensors was fixed (in both papers, three sensors are used). The aim of an analysis of the numerical model (1) is to determine position(s) of sensor(s) to identify thermal properties of the material as accurately as possible.
D. Jarušková, A. Kucˇerová / International Journal of Heat and Mass Transfer 106 (2017) 135–141
136
Nomenclature thermal conductivity along the xaxis thermal conductivity along the yaxis true value of kx true value of ky thermal volumetric capacity temperature field sample dimensions u heat flux density ru standard deviation of u J number of sensors x x-coordinate of a sensor position y y-coordinate of a sensor position DxðjÞ; j ¼ 1; . . . ; J random displacements along the x-axis DyðjÞ; j ¼ 1; . . . ; J random displacements along the y-axis rx standard deviation of DxðjÞ ry standard deviation of DyðjÞ n number of measurements t i ; i ¼ 1; . . . ; n time points of measurement
kx ky kx ky C T Lx ; Ly
Ruffio et al. [4] claim that a solution of (1) may be found analytically in the form of a quickly converging series. In our paper, we assume that a final solution may be quite accurately approximated by the first term only, i.e.,
Tðx; y; tÞ ¼ hx ðt; kx ; C; u; xÞ þ hy ðt; ky ; C; u; yÞ;
where
u pffiffi p hx ðt; kx ; C; u; xÞ ¼ p2ffiffiffiffiffi t F exffit ; Ckx e u ffi pffiffi p t F yffi hy ðt; ky ; C; u; yÞ ¼ p2ffiffiffiffiffi t
Cky
pffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffi with e x ¼ ðx=2Þ C=kx , e y ¼ ðy=2Þ C=ky and
FðzÞ ¼
Z expðz2 Þ 2 pffiffiffiffi z 1 pffiffiffiffi
p
p
0
z
2 ev dv ;
z P 0:
Because sensors measure temperature with measurement errors, we may assume that the behavior of a measured temperature can be modelled using a nonlinear regression and the problem described above is a design experiment problem in the nonlinear regression [6]. In the case of three sensors, it was numerically solved in [5]. However, Ruffio et al. [4] argue that measured temperature is not only affected by measurement errors but other input parameters may be random as well. (The very same problem was also considered by Emery et al. [1]). For instance, heat flux density applied to sample boundaries is not perfectly known. Similarly, a sensor may not be located exactly at a designed position. Then, instead of being measured at a point ðx; yÞ temperature is measured at a point ðx þ Dx; y þ DyÞ, where Dx and Dy are random. Ruffio et al. [4] assume that the distribution of the variables u; DxðjÞ; DyðjÞ; j ¼ 1; . . . ; J (J being a number of sensors), is known. In this way, they deal with the design of an experiment in a nonlinear regression with random parameters. Additionally, it is implicitly assumed that random parameters do not change during one experiment. In the other words, each of them attains only one fixed value for all time points t ¼ 1; . . . ; 60, but they may attain different values if the experiment is repeated. This means that these random parameters play the role of so-called random common factors. Then, as we illustrate below, the variability of the estimates of the parameters of interest is practically entirely determined by the variability of random parameters, i.e. it depends on the variability of measurement errors only very slightly. The goal of an analysis is to find the ‘‘optimal” position(s) of sensor(s) so that if
s final time of measurement Y i ; i ¼ 1; . . . ; n values of the dependent variable b vector of the parameters of interest b true value of b c vector of random parameters c0 mean of c Rc variance–covariance matrix of c n design parameter rm standard deviation of measurement errors f ; g; h regression functions F matrix of partial derivatives of f with respect to the components of b computed at values b and c0 Z matrix of partial derivatives of f with respect to the components of c computed at values b and c0 v21a ðpÞ ð1 aÞ 100% quantile of the v2 distribution with p degrees of freedom MSEðbÞ mean square error
we repeat the experiment, the differences between estimates of parameters and their true values are small. If we planned to do only one experiment, it would be better to consider the random parameters to be nuisance parameters that might be estimated. Their known distribution might serve as a prior distribution and we might try to find the estimates of the parameters of interest as well as of the nuisance parameters using a Bayesian approach. This is, however, a different statistical problem than we are solving here. Despite general agreement with the method for designing an experiment as suggested by Ruffio et al. [4], we are afraid that the authors do not fully realize problems connected to nonlinear regression models with random factors. Our objections do not only concern the specific problem described above, but all similar problems. For the sake of clarity, we simplified the situation assuming that the parameters of interest are heat conductivities kx and ky only, while thermal volume capacity C is known to be equal to 1,700,000 (Jm3 K1 ). 3. Non-linear regression with non-random parameters We suppose that we observe a sequence of pairs fðt i ; Y i Þ; i ¼ 1; . . . ; ng. We assume that ft i ¼ is=n; i ¼ 1; . . . ; ng are time points from the interval ½0; s and
Y i ¼ f ðt i ; bÞ þ ei ;
i ¼ 1; . . . ; n;
ð2Þ
where f ðt; bÞ is a known regression function of t 2 ½0; s and a p-dimensional vector of unknown parameters of interest b ¼ ðb1 ; . . . ; bp ÞT . The random errors fei g are i.i.d. random variables with a variance r2m . We suppose that the true value of b is b 2 ½bl ; bu and we estimate b using the least squares method: n X 2 b b ¼ argmin ðY i f ðt i ; bÞÞ : b2½bl ;bu i¼1
Under a set of general conditions, b b is a consistent estimate of b pffiffiffi P (i.e. b b ! b as n ! 1) and nð b b b Þ has asymptotically a 1 p dimensional normal distribution N 0; ðF T F Þ r2m , see [6]. The matrix F is a matrix of partial derivatives of f with respect to b at the right value b .
D. Jarušková, A. Kucˇerová / International Journal of Heat and Mass Transfer 106 (2017) 135–141
Further, we assume that the regression function f ðt; b; nÞ depends on some parameter n. The goal of our analysis is to choose n in such a way that the convergence of b b to b is as fast as possible. As the volume of a ð1 aÞ100% asymptotic confidence ellipsoid
T ðb b b Þ F T F ð b b b Þ 6 r2m v21a ðpÞ is proportional to
of b b under an assumption that at the same time the parameter c has been randomly chosen from a distribution N c0 ; r2c and the measurement errors fei ; i ¼ 1; . . . ; ng have been randomly chosen b by from Nð0; r2 Þ. Here, we can express b m
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 det F T F . The D criterion consists in mini-
mizing a determinant of an asymptotic variance–covariance matrix 1 . We will of b b and that is here the same as to minimize det F T F also consider the F criterion, introduced by Ruffio et al. [4], where P 1=2 p 2 b is to be minimized. i¼1 Var b i =ðbi Þ 1 The asymptotic variance–covariance matrix F T F r2m depends on the right value b so that an optimal value nop with respect to the D as well as the F criterion depends on b as well. This is a substantially different than in linear regression case where a choice of nop is independent of b , see [3]. One way to overcome the problem is to proceed in two steps, i.e. in the first step, estimating the parameter b only roughly and in the second step, finding nop with respect to the estimate obtained in the first step. If we know that b 2 ½bl ; bu , we can minimize a criterion with respect to the ‘‘worse” b in the interval ½bl ; bu . This is basically the idea of Ruffio et al. [4]. 4. Non-linear regression with random parameters
P is is P is h n ;n g n ;n g ;n e b P 2 is b ¼ b þ c c0 þ P 2nis i : g n ;n g n ;n
As c and fei g are independent, the unconditional distribution of b b is normal with the mean b and a variance:
P is is !2 h n ;n g n ;n r2 P 2 is r2c þ P 2 mis : g n ;n g n ;n
Var b b¼
i ¼ 1; . . . ; n: ð3Þ
: b¼ Var b
Rs 0
hðt; nÞgðt; nÞdt Rs g 2 ðt; nÞdt 0
!2
r2c þ
r2m n
Rs 0
1 : g 2 ðt; nÞdt
may minimize either (4) or (5). For a large n, the difference will be small. In the case of a nonlinear regression model with a random parameter c:
i ¼ 1; . . . ; n;
ð6Þ
Ruffio et al. [4] suggest replacing the model (6) with a linear model
Y i ¼ f ðt i ; b ; c0 ; nÞ þ þ
@f ðt i ; b ; c0 ; nÞðb b Þ @b
@f ðt i ; b ; c0 ; nÞðc c0 Þ þ ei ; @c
i ¼ 1; . . . ; n
@f ðt ; b ; @c i
parameter c is independent of random errors fei ; i ¼ 1; . . . ; ng that are i.i.d. random variables. For simplicity, we suppose that they are distributed according to Nð0; r2m Þ. Suppose that in our experiment, the parameter c actually c that might be different from c0 . As we do not attained a value e Pn know this value, we estimate b by minimizing i¼1 ðY i l
1 1 1 F T F F T Z Rc Z T F F T F þ r2m F T F ;
b b ¼ b þ
Pn
is i¼1 ðY i Þg n ; n Pn 2 is ;n i¼1 g n
l
:
The conditional distribution of b b given c ¼ e c describes a stochastic behavior of b b under the assumption that the experiment is perc . It is normal with formed many times under a condition that c ¼ e a mean
b þ
P is is h n ;n g n ;n e P 2 is c c0 g n ;n
@f This is indeed model (3) with gðt i ; nÞ ¼ @b ðt i ; b ; c0 ; nÞ and hðt i ; nÞ ¼
c0 ; nÞ.
For a multidimensional parameter b, we obtain an approximate b in the form presented by Ruffio variance–covariance matrix of b et al. [4]:
r2m
g 2 ins ; n
:
The unconditional distribution of b b under an assumption that c has a normal distribution N c0 ; r2c describes a stochastic behavior
ð7Þ
where F is a matrix of partial derivatives of f with respect to components of b, Z is a matrix of partial derivatives of f with respect to components of c, both computed at b and c0 , and Rc is a diagonal matrix with a diagonal ðVarc1 ; . . . ; Varck Þ. Ruffio et al. [4] suggest applying either the D criterion or the F criterion to the matrix (7). Similarly as in the linear regression with random parameters, b is a non-consistent estimator. Moreover, b b is no the estimate b longer an unbiased estimator and its bias does not disappear as the number of observations tends to infinity. The reason is the b of the least squares following. We are looking for a minimum b function MSEðbÞ:
MSEðbÞ ¼
and a variance
P
ð5Þ
The estimate b b is an unbiased estimator but it is not consistent Rs unless 0 hðt; nÞgðt; nÞdt ¼ 0. The asymptotic distribution of b b b is a zero mean normal distribution with a variance equal to the first term in (5). In case we wish to minimize Var b b with respect to n, we
We suppose that gðt; nÞ and hðt; nÞ are continuous functions on ½0; s ½nl ; nu . We suppose that l is a known real parameter, while b is a parameter that is to be estimated. Further, we suppose that c is a random parameter distributed according to Nðc0 ; r2c Þ. The
gðis=n; nÞðb b ÞÞ2 . Therefore, the sought estimate has a form:
ð4Þ
If n is large, the variance Var b b behaves approximately as
Y i ¼ f ðt i ; b; c; nÞ þ ei ; In this section we explain properties of a least squares estimate of a parameter of interest in a non-linear regression with random parameters. We start analyzing a linear regression model in the form:
Y i ¼ l þ gðis=n; nÞðb b Þ þ hðis=n; nÞðc c0 Þ þ ei ;
137
¼
n 2 1X Y i f ðt i ; b; c0 ; nÞ n i¼1 n 2 1X f ðti ; b ; c; nÞ f ðt i ; b; c0 ; nÞ n i¼1
þ
n n 2X 1X ei ðf ðt i ; b ; c; nÞÞ f ðt i ; b; c0 ; nÞÞ þ e2 : n i¼1 n i¼1 i
ð8Þ
D. Jarušková, A. Kucˇerová / International Journal of Heat and Mass Transfer 106 (2017) 135–141
138
As n tends to infinity, the second term of (8) converges to zero and Rs b minimizes ðf ðt; b ; c; nÞ the third term to r2 . The optimal b m
0
f ðt; b; c ; nÞÞ dt and it is generally a nonlinear function of c. Due to b may differ from b , nonlinearity, the mean of the distribution of b 0
2
b may differ from the matrix and the variance–covariance matrix of b (7). Clearly, this discrepancy does not disappear as n tends to infinb canity. Therefore, an approximate variance–covariance matrix of b b unless not be a criterion for determining accuracy of the estimate b it is approximately a linear function of c. For a one dimensional parameter, we should compare the quality of an estimate b b using 2 a square mean deviation from the true value E ð b b b Þ . We recall 2 2 b þ ðE b that Eð b b b Þ ¼ Var b b b Þ . For a p-dimensional paramePp 2 2 b ter, we consider the F criterion, where is to j¼1 Eð bj bj Þ =bj
be minimized. However, we cannot evaluate it unless we know the true value b and we get into a vicious circle.
we infer the matrix F comprising of two columns, namely x ¼ ðx1 ; . . . ; xn ÞT and y ¼ ðy1 ; . . . ; yn ÞT , where for i ¼ 1; . . . ; n
pffiffiffiffi e 2u x xi ¼ pffiffiffi 3=2 ti G pffiffiffiffi ; ti C ðkx Þ with e x ¼ ðx=2Þ
pffiffiffiffi e 2u y yi ¼ pffiffiffi 3=2 t i G pffiffiffiffi ti C ðky Þ
qffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffi y ¼ ðy=2Þ C=ky and kx , ky being true values of C=kx , e
parameters kx , resp. ky . First, we consider the D criterion. Without a loss of generality we confine ourselves to a special case of kx ¼ ky ¼ 1 and search b y op ð1; 1Þ as for general values of the conductivities kx x op ð1; 1Þ; b and ky . It holds that
b x op ðkx ; ky Þ ¼
qffiffiffiffiffi kx b x op ð1; 1Þ;
b y op ðkx ; ky Þ ¼
qffiffiffiffiffi ky b y op ð1; 1Þ
ð10Þ
and for
1
T ; Dop ð1; 1Þ ¼ min det ðF ð1; 1ÞÞ ðF ð1; 1ÞÞ x;y
5. Numerical results
(
In this section we highlight some interesting properties of the design of an experiment in a nonlinear regression. We study model (1) with Lx ¼ Ly ¼ 0:05 (m), see Fig. 1. Recall that the unknown parameter of interest is the vector of thermal conductivities ðkx ; ky ÞT . We start with a non-linear regression where no random factor is present. Case 1 We are to find a position ðx; yÞ for one sensor to estimate b ¼ ðkx ; ky ÞT . We assume that all other parameters are fixed; the
Dop
¼ min det
x;y
F
ðkx ; ky Þ
T
F
ðkx ; ky Þ
) 1 ;
it holds
3 Dop ¼ kx ky Dop ð1; 1Þ:
ð11Þ
We proved analytically that for kx ¼ ky ¼ 1 there exist two optimal symmetric positions that are situated on the boundaries of the sample. Looking for an optimal value of the x-coordinate of the first optimal position on a grid ½0 : 0:0001 : 0:05 we have found that ð1Þ
ð1Þ
volumetric thermal capacity C ¼ 1; 700; 000 ðJm3 K1 Þ, the heat
ðxop ¼ 0:0023; yop ¼ 0Þ. Clearly, for the second optimal position it
flux density u ¼ 25; 000 ðWm2 Þ, and the sensor is located exactly at ðx; yÞ so that Dx ¼ 0 and Dy ¼ 0. In agreement with [4], temperature is measured at times t i ¼ i; i ¼ 1; . . . ; n ¼ 60. A standard deviation of measurements errors rm ¼ 0:1. Under these assumptions, the regression function’s form is:
holds
e e 2u pffiffi 2u pffiffi x y f ðt; kx ; ky ; x; yÞ ¼ pffiffiffiffiffiffiffiffi tF pffiffi þ pffiffiffiffiffiffiffiffi t F pffiffi : Ckx Cky t t
ð9Þ
For finding an optimal ðx; yÞ, we compute the asymptotic variance– covariance matrix of b b. Introducing a new function GðzÞ by
Z 2 GðzÞ ¼ FðzÞ þ z 1 pffiffiffiffi
p
0
z
2 ev dv ;
ð2Þ ðxop
ð2Þ
¼ 0; yop ¼ 0:0023Þ. The minimal value of the D criterion 13
is 1:40 10 . For kx ¼ 0:6 and ky ¼ 4:7, considered in [4], the optimal positions of a sensor and the optimal value of the D criterion are pffiffiffiffiffiffiffi ð1Þ ð1Þ ð2Þ ð2Þ indeed xop ¼ 0:6 0:0023 ¼ 0:0018; yop ¼ 0 and xop ¼ 0; yop ¼ pffiffiffiffiffiffiffi 4:7 0:0023 ¼ 0:0050, Dop ¼ ð0:6 4:7Þ3 1:40 1013 ¼ 3:14 1012 . Notice, see (10), that an optimal position may be substantially affected by the true values of parameters of interest. The criterion F is not homogeneous with respect to kx and ky and that is why we have to find the optimal positions for kx ¼ ky ¼ 1 and kx ¼ 0:6; ky ¼ 4:7 separately. For kx ¼ ky ¼ 1 the optimal posi ð1Þ ð1Þ ð2Þ ð2Þ tions are xop ¼ 0:0025; yop ¼ 0 and xop ¼ 0; yop ¼ 0:0025 with an optimal value 1:8 103 . We see that in this case and kx ¼ ky ¼ 1 the optimal positions for the criteria D and F are not the same but they differ only slightly. For kx ¼ 0:6 and ky ¼ 4:7 only one optimal solution xop ¼ 0:0024; yop ¼ 0 exists and the optimal value of the criterion F is 0:0016. For kx ¼ ky ¼ 1 the matrix ðF ð1; 1ÞÞ ðF ð1; 1ÞÞ is singular when a sensor is placed on the diagonal of the square sample, i.e. x ¼ y, as the both coordinates of F ð1; 1Þ are the same. Positions near the diagonal yield an ill-conditioned matrices, e.g. the ratio of eigenT
ð1Þ
ð1Þ
values for ðxop ¼ 0:0023; yop ¼ 0Þ is 75.8. For a general kx ; ky the qffiffiffiffiffi pffiffiffiffiffi matrix F T F is singular when x= kx ¼ y= ky . To avoid the problem with singularities instead of minimizing the D and F criterion we maximized their inverses. In the three following cases, we estimate the parameter
Fig. 1. Description of experiment to optimize.
b ¼ ðkx ; ky ÞT in the nonlinear regression model (6) with the regression function (9) and a standard deviation of measurement errors rm ¼ 0:1, but we suppose that some of the input parameters are random. We compare the value of the F criterion for the ‘‘optimal” sensor position where the variances were obtained using the lin-
D. Jarušková, A. Kucˇerová / International Journal of Heat and Mass Transfer 106 (2017) 135–141
earization with the value of the F s , resp. F s , criterion where means and variances were obtained by numerical solutions of a non-linear minimization problem with parameters obtained by the Monte Carlo simulations. Our aim is to show that in some cases ‘‘good design” obtained by linearization yields very bad design as the mean square errors of the least squares estimates of the identified parameters are large. The optimal position(s) of sensor(s) were found numerically on a grid ½0 : 0:0001 : 0:05 ½0 : 0:0001 : 0:05 or using an evolutionary GRADE algorithm extended by niching strategy CERAF [2]. The matrices used in the D and F criteria were the approximate variance–covariance matrices given in (7). After finding the optimal sensor position, we simulated 100,000 nonlinear models (6) with random common factors. Using the least squares method, we estimated the parameter of interest b for all 100,000 simulated models. For obtaining the least squares estib we applied the Matlab function ‘‘lsqcurvefit” with minimates b, mum precision of the target least squares function ‘TolFun’ equal to 1025 and minimum precision of the estimated parameters ‘TolX’ equal to 1015 . Case 2 We consider the nonlinear regression model (6) with the regression function (9) and the true values of parameters kx ¼ ky ¼ 1. We suppose that heat flux density u is a random variable distributed according to a normal distribution with a mean
u0 ¼ 25; 000 and ru ¼ 50 ðWm2 Þ. All other parameters are fixed.
We are looking for the position of one sensor ðx; yÞ. First, we consider the D criterion. Fig. 2 shows contours of detðVar b bÞ on ½0; 0:005 ½0; 0:005 computed by (7) with respect
to ðx; yÞ. The optimal value is attained at two positions x ¼ 0:0027; y ¼ 0 and x ¼ 0; y ¼ 0:0027 and the attained minimum is 3:04 1011 . Using the results of simulations for x ¼ 0:0027; y ¼ 0 k y 1Þ ¼ 0:0006 and we estimate the bias Eðb k x 1Þ ¼ 0:0002; Eðb k y ¼ 0:0008. the variances Var b k x ¼ 0:0002; Var b For the F criterion, the optimal value is attained at two positions x ¼ 0:0023; y ¼ 0 and x ¼ 0; y ¼ 0:0023. The attained minimum is 0:0036. Using the results of simulations for x ¼ 0:0023; y ¼ 0, we estimate the bias Eðb k x 1Þ ¼ 0:0002; Eðb k y 1Þ ¼ 0:0006, and the k y ¼ 0:0008 and the correlation coefvariances Var b k x ¼ 0:0002; Var b b k y Þ ¼ 0:96. ficient Corrð k x ; b We have seen that the accuracy of the parameters’ estimates in regression with random parameters is given by the variability of the random factors. Because the variability of the heat flux density is small, the optimal sensor positions for a random and nonrandom u are very close. Moreover, the dependence of b k x and kby on u is practically linear and therefore the biases of the estimates are negligible. This is illustrated by Fig. 3 that shows the dependence of b k x k and kby k on u 2 ½u0 3ru ; u0 þ 3ru for x
y
x ¼ 0:0023; y ¼ 0; kx ¼ 1 and ky ¼ 1. −3
5
y.10
Fig. 3. Dependence of the estimates b k x kx (solid line) and b k y ky (dashed line) on ufor kx ¼ ky ¼ 1 and x ¼ 0:0023; y ¼ 0.
Case 3 We consider again the nonlinear regression model (6) with the regression function (9) with the true values of parameters kx ¼ 0:6 and ky ¼ 4:7. We would like to measure temperature using one sensor and we know that we are not able to place it exactly at the designed position. The random displacements Dx and Dy from the right position ðx; yÞ are i.i.d., distributed according to a normal distribution with zero mean and a standard deviation rx ¼ ry ¼ 0:0005 (m). All other parameters are fixed. Using the F criterion, the optimal solution is attained at x ¼ 0:0165; y ¼ 0:0158 and the minimum of the criterion F is 0.0387. The contours of the F criterion are presented in Fig. 4. Notice that the ‘‘optimal” position differs substantially from the optimal position when the random common factors are not present. Using simulated values, the biases estimates are Eðb k x 0:6Þ ¼ 0:0001 and Eðb k y 4:7Þ ¼ 0:0436. The variances are k y ¼ 0:2486 and Corrðb kx ; b k y Þ ¼ 0:87. The value Var b k x ¼ 0:0117; Var b of the F s criterion is equal to 0.2093 and that of F s to 0.2095. We see that the biases are small but the values of F s , resp. of F s criterion, obtained by simulations, differ from the values of the F criterion obtained using (7). Finally, we simulated model (6) without measurement errors, i.e. r2m ¼ 0. For x ¼ 0:0165 and y ¼ 0:0158, the estimates of the biases and the variances changed only very little. The biases are 0.0000 and 0.0403; the variances are 0.0111 and 0.2483. We see that if n is large, the random measurement errors do not affect variability of the parameters’ estimates. Case 4 In this case we compare the quality of four designs suggested by Ruffio et al. [4] and denoted by U1; U2; U3, and U4. All four designs correspond to the location of three sensors. The aim of the analysis presented in [4] was to find the optimal positions not only to estimate thermal conductivities but also the volumetric thermal capacity C. In accordance with the previous cases, we assess the designs according to a quality of estimates of thermal conductivities kx and ky assuming C ¼ 1; 700; 000 Jm3 K1 . Further, in agreement with [4] we assume that u is a normally distributed random y [m] 0.03
[m]
4
0.025
3
0.02
2
0.015
1 0
139
0
1
2
3
4 5 −3 x.10 [m]
Fig. 2. Contours of det Var b b computed using an approximate variance–covariance matrix (7) for ðx; yÞ 2 ½0; 5 103 ½0; 5 103 .
0.01 0.01
0.015
0.02
0.025
0.03 x [m]
Fig. 4. Contours of criterion F computed using an approximate variance–covariance matrix (7) for ðx; yÞ 2 ½0:01; 0:03 ½0:01; 0:03.
D. Jarušková, A. Kucˇerová / International Journal of Heat and Mass Transfer 106 (2017) 135–141
140
variable with a mean u ¼ 25; 000 and a standard deviation
ru ¼ 50 (W m2 Þ, and that, for j ¼ 1; . . . ; 3, the variables DxðjÞ as
well as DyðjÞ had zero mean normal distribution with a standard deviation of 0.0005. All random variables are independent. The true value of kx is 0.6 and of ky is 4.7. Paper [4] does not provide the exact positions of the sensors for U1; U2; U3, and U4 and we interfered them from Fig. 10, [4]. We assume the positions to be: U1
U3
x1 ¼ 0:011 x2 ¼ 0:011 x3 ¼ 0:037
y1 ¼ 0:031 y2 ¼ 0:031 y3 ¼ 0:05
U2
x1 ¼ 0:010 x2 ¼ 0:033 x3 ¼ 0:044
y1 ¼ 0:031 y2 ¼ 0 y3 ¼ 0
U4
x1 ¼ 0:033 x2 ¼ 0:044 x3 ¼ 0:016
y1 ¼ 0:023 y2 ¼ 0:023 y3 ¼ 0:05
x1 ¼ 0:011 x2 ¼ 0:030 x3 ¼ 0:039
y1 ¼ 0:048 y2 ¼ 0 y3 ¼ 0
We can see that all designs have either one or two sensors on the lower or upper boundary of a sample. In such a case, the corresponding Dy cannot have a zero mean normal distribution. Therefore, we assume that if a sensor was located on the lower boundary y ¼ 0, then Dy had a (positive) half normal distribution with a standard deviation stdðDyÞ ¼ 0:0005. If it was located on the upper boundary, then Dy had a (negative) half normal distribution with a standard deviation stdðDyÞ ¼ 0:0005. In such a way, the mean of ycoordinate of a sensor originally located on the lower boundary is 6:6 104 , and of a sensor originally located on the upper boundary is 0:05 6:6 104 . Table 1 presents values of F based on the approximate variances obtained by (7) and the values of F s and F s obtained by simulations. All biases are small. However, by ordering the designs according to the performance of F based on linearization, we get: U2; U1; U4, and U3 while by ordering according to the performance of F s , resp. F s , we get: U2; U4; U1; U3. This is caused by the difference between the values of F and F s . Using a genetic GRADE algorithm extended by niching strategy CERAF [2] we found an ”optimal” design U5 with respect to the F criterion:
x1 ¼ 0:025327;
y1 ¼ 0:024367;
x2 ¼ 0:002705;
y2 ¼ 0:012212;
x3 ¼ 0:004806;
y3 ¼ 0:005897:
Using (7) the value of the F criterion is equal to 0.06 as Var b k x ¼ 0:0006 and Var b k y ¼ 0:0449. Using simulations we estib k y equal to mated the bias of k x equal to 0.129 and the bias of b k y ¼ 13:13; F s ¼ 0:973 1.078. Moreover, we get Var b k x ¼ 0:127; Var b
and F s ¼ 1:019. Clearly, the design U5 is a bad design comparing it with U1; U2; U3, and U4 despite the value of the F criterion computed by (7) being very small. Not only the estimates have extremely large biases but the variances are large as well. Fig. 5 illustrates nonlinearity of relationship between b k x , resp. b k y , and
the x-coordinate of the third sensor. We would like to remark that for all analyzed designs a difference between the corresponding elements of the variance–
Table 1 The values of F based on the approximate variances obtained by (7) and the values of F s and F s obtained by simulations. U1 U2 U3 U4
F F F F
¼ 0:1285 ¼ 0:0964 ¼ 0:1649 ¼ 0:1478
Fs Fs Fs Fs
¼ 0:1470 ¼ 0:0877 ¼ 0:1549 ¼ 0:1400
F s F s F s F s
¼ 0:1473 ¼ 0:0883 ¼ 0:1546 ¼ 0:1404
Fig. 5. Dependence of the estimates ðb k x kx Þ=kx (solid line) and ðb k y ky Þ=ky (dashed line) on x3 for kx ¼ 0:6; ky ¼ 0:8 and the design U5.
covariance matrix for the model with and without measurement errors is less 103 . Finally, we study properties of the robust version of the criterion F. If it is known that the parameters of interest are elements of some known intervals, i.e.,
b 2 ½b1 db1 ; b1 þ db1 . . . ½bp dbp ; bp þ dbp ; the authors [4] suggest evaluating the F criterion not only for b ¼ ðb1 ; . . . ; bp ÞT but also for fðb1 db1 ; . . . ; bp dbp ÞT g and compute the maximum of all values of F as a worst case scenario. We supposed that kx 2 ½0:3; 0:7 and ky 2 ½3:0; 7:0 and computed the values of F s and F s obtained by simulations for all designs U1; U2; U3, and U4 and ðkx ; ky Þ ¼ ð0:6; 4:7Þ; ð0:3; 3:0Þ; ð0:3; 7:0Þ; ð0:7; 3:0Þ; ð0:7; 7:0Þ. The designs U1; U3 and U4 are relatively unsensitive to the change of true values of kx and ky , but the design U2 is sensitive to the change of kx . Applying a genetic GRADE algorithm [2] we have found the ‘‘optimal design” U6 with the respect to FR criterion:
x1 ¼ 0:040665;
y1 ¼ 0:034156;
x2 ¼ 0:014721;
y2 ¼ 0:046086;
x3 ¼ 0:014148;
y3 ¼ 0:05:
The value of F is 0.089, and of FR is 0.091. Using simulations we found that the biases for kx 2 ½0:3; 0:7 and ky 2 ½3:0; 7:0 are generally small. The value of F s is 0.073 and of its robust version is 0.160. Therefore, we consider the U6 design to be very good. 6. Conclusion We numerically analyzed a heat transfer problem described by Sawaf et al. [5] and Ruffio et al. [4] with the aim of finding optimal sensor positions for estimating thermal conductivities. We have shown that if we use the approach suggested in [4] with a criterion based on an approximate variance–covariance matrix of estimates in nonlinear regression models with random parameters, we encounter some problems. First, if the number of measurements n is large, the level of inaccuracy of estimated parameters is almost completely determined by the variability of random parameters (i.e. common factors). Moreover, the estimates of the parameters of interest are not consistent, i.e. the estimate does not tend to the true value as the number of measurements increases. In such a case, their possible bias that arises from a nonlinear dependence of the estimates of parameters of interest on random ‘‘nuisance” parameters may play an important role. Non-linearity may also cause the substantial difference between the variance–covariance matrix obtained using linearization from the true variance–covariance matrix of the estimates. In some problems, approximate linearity may be analyzed by repeating numerical calculations for those parameters
D. Jarušková, A. Kucˇerová / International Journal of Heat and Mass Transfer 106 (2017) 135–141
values that are attained most ‘‘likely”. An alternative method is to assess the quality of estimates using Monte Carlo simulations. It would be difficult to find an optimal design with simulations, but at least the characteristics of the solution obtained using linearization should be compared with the characteristics obtained via simulations. References [1] A.F. Emery, A.V. Nenarokomov, T.D. Fadale, Uncertainties in parameter estimation: the optimal experiment design, Int. Heat Mass Transfer 43 (18) (2000) 3331–3339.
141
[2] A. Ibrahimbegovic, C. Knopf-Lenoir, A. Kucˇerová, P. Villon, Optimal design and optimal control of structures undergoing finite rotations and elastic deformations, Int. J. Numer. Meth. Eng. 61 (2004) 2428–2460. [3] L. Pronzato, A. Pázman, Design of Experiments in Nonlinear Models, Springer, New York, 2013. [4] E. Ruffio, D. Saury, D. Petit, Robust experiment design for the estimation of thermophysical parameters using stochastic algorithms, Int. Heat Mass Transfer 55 (2012) 2901–2915. [5] B. Sawaf, M.N. Ozisik, Y. Jarny, An inverse analysis to estimate linearly temperature dependent thermal conductivity components and heat capacity of an orthotropic medium, Int. Heat Mass Transfer 38 (16) (1995) 3005–3010. [6] G.A.F. Seber, C.J. Wild, Nonlinear Regression, Wiley, New York, 1989. [7] A. Tarantola, Inverse Problem Theory and Methods for Model Parameter Estimation, Siam, Philadelphia, 2005.