A D-Optimal Sequential Calibration Design for Computer Models

: The problem with computer model calibration by tuning the parameters associated with computer models is signiﬁcant in many engineering and scientiﬁc applications. Although several methods have been established to estimate the calibration parameters, research focusing on the design of calibration parameters remains limited. Therefore, this paper proposes a sequential computer experiment design based on the D-optimal criterion, which can efﬁciently tune the calibration parameters while improving the prediction ability of the calibrated computer model. Numerical comparisons of the simulated and real data demonstrate the efﬁciency of the proposed technique.


Introduction
Experiments are conducted to explore or optimize physical phenomena. In some applications, such as national defense, medicine, and manufacturing, the physical experiments may be difficult to be conducted due to economic, technical, or ethical limitations. To reduce the experimental cost, mathematical models, which are also called computer models, are developed to mimic, understand, and predict the physical phenomena in many applications [1][2][3]. The computer models are useful and efficient only if they can approximate the physical process well. Oftentimes, the computer models contain a set of calibration parameters, which are physical unobservable variables. The computer models' fidelity to physical process relies on the unknown values of calibration parameters. Then, physical data and computer model outputs are combined to estimate the calibration parameters such that the computer model matches the physical process. This procedure is referred to as computer model calibration in the literature.
Numerous models have been proposed in the literature for the problem of computer model calibration, such as [2,[4][5][6][7]. Among them, the Kennedy and O'Hagan model is the most commonly used. The Kennedy and O'Hagan model integrates the physical data and computer model outputs through a Bayesian framework. Any posterior quantity can serve as the point estimate of calibration parameter depending on the loss function specified. In practice, the most commonly used predictor of the physical process is a calibrated computer model, see [8,9]. Since the nonlinear effects from the discrepancy function can be hard to interpret and also may open up the possibility of overfitting with limited physical observations. Ref. [10] pointed out that an interpretable calibration parameter should allow the computer model to predict the real physical phenomena well even without the discrepancy function.
How to perform the experiments efficiently to tune the calibration parameters accurately under some metrics plays an important role. Although we are not the first to look into the problem of design for calibration, it has not received enough attention. Based on the Kennedy and O'Hagan model, some designs have been proposed in the literature. Ref. [11] employed the Kullback-Liebler (KL) divergence criterion as a function of the computer model inputs and obtained the estimate of calibration parameters by minimizing this criterion. Ref. [12] focused on the problem of functional calibration by generating sequential designs for the physical and computer experiments. Ref. [13] used results from the nonlinear optimal design theory to design such experiments. Ref. [14], based on the Kennedy and O'Hagan model, proposed an optimal sequential design for both computer and physical experiments by regarding integrated mean squared prediction error. Based on the Bayesian model calibration framework of [6], a D-optimal design for the physical experiment was proposed by [15]. Ref. [16] proposed a follow-up optimal experiment design for computer models calibration. In some practices, no physical experiments can be conducted after the initial design due to limitations. Thus, these designs, considering the physical experiments, are hard and impracticable. Ref. [17] presented an adaptive design for computer experiments to estimate the calibration parameters by using the expected improvement (EI) algorithm. It aims at reducing the calibration error induced by the uncertainty of the emulator of computer models but not at improving the estimation of calibration parameters. Inspired by this, we divert effort on the designs only considering computer experiments to estimate the calibration parameters. The D-optimal criterion is well known and widely used in the literature, which can help gather more information about the calibration parameters by minimizing the asymptotic variances of estimate. This paper proposes a sequential computer experiment calibration design using the D-optimal criterion and presents a fast algorithm to generate the designs.
The article is organized as follows. Section 2 reviews the Kennedy and O'Hagan calibration method. In Section 3, the proposed local D-optimal sequential design is presented. A fast algorithm for generating the corresponding designs is suggested. In Section 4, some simulation studies are made to demonstrate the performance of the proposed design. Conclusions and remarks are given in Section 5. Appendix A shows the derivation of the Fisher Information Matrix (FIM) for the calibration parameters.

Calibration of Computer Models
An important reference for computer calibration is the work of [6]. In this section, we will review some related background about the Kennedy and O'Hagan model. Let y p be the observation of physical process and x = (x 1 , · · · , x r ) ∈ X ⊂ R r be the control variables, which are also the set of inputs for physical process. According to [6], the physical observation y p can be modeled as where y c (·, ·) is a computer model, θ ∈ T ⊂ R h is the set of calibration parameters, δ(·) is the discrepancy function which is independent of the computer model y c (·, ·); ∼ N(0, λ 2 ) is the observation error and λ 2 is the corresponding variance. In the literature, the most popular methods to fit the computer model y c (·, ·) and discrepancy function δ(·) are the Gaussian processes due to analytical tractability. Thus, the prior information about both y c (·, ·) and δ(·) is considered as y c (·, ·) ∼ GP(m(·, ·), c 1 {(·, ·), (·, ·)}), where m(·, ·) is the mean function of computer model. Assume x and x denote the values of control inputs, and t and t denote the values of calibration inputs. According to the literature [18,19], c 1 {(·, ·), (·, ·)} and c 2 (·, ·) are usually the corresponding separable covariance functions such as Here, σ 2 1 and σ 2 2 are the variance parameters, and In terms of the mean function m(·, ·), the linear model structure is always considered, i.e., where h(x, t) = (h 1 (x, t), h 2 (x, t), · · · , h p (x, t)) T is a vector of p known functions over X and β = (β 1 , β 2 , · · · , β p ) T is the corresponding unknown regression coefficients.
p q } be the design for physical experiment with q points, y p = {y p 1 , y p 2 , · · · , y p q } T be the corresponding physical outputs, be the design for computer experiment with n points, and y c = {y c 1 , y c 2 , · · · , y c n } T be the corresponding computer outputs. Thus, the full output d = (y cT , y pT ) T is normally distributed given (θ, β, σ 2 1 , σ 2 2 , λ 2 , Ω x , Ω θ , Ω δ ), and the corresponding likelihood function can be yielded. In order to express the mean and variance matrix of full output clearly, we define the following notations. Let ψ 1 = (σ 2 1 , , θ)} be the augmented design points by calibration parameters θ. Then the mean and variance matrix for full output vector d given (θ, β, ϕ) can be derived as and where ; and I q is a q × q identity matrix. Then the posterior for the parameters given the data can be written as where π(θ, β, ϕ) is the prior for unknown parameters. The MCMC techniques are usually used to determine the posterior distribution, but require complex computations. To simplify the computations, we adopt the modularization by the literature [20], namely, first estimate the emulator of the computer model and then the discrepancy.
The modular approach in [20] is considered here to estimate the parameters in the model, which can be described as follows. The maximum likelihood estimates (MLEs)β of β andψ 1 of ψ 1 can be obtained based on the computer experimental data (D c , y c ). For the calibration parameter θ, which is a tuning parameter, there is no "true value". The goal of calibration is to find out some type of best-fitting value of θ. It is easy to obtain the least-squares estimate of θ, i.e., the valueθ that minimizes the differences between the physical outputs and the computer model outputs by regarding fixing the β and ψ 1 at their MLEs. The bias data (D p , y p −ŷ c (D p θ )) can be employed to compute the MLEsψ 2 of ψ 2 andλ of λ, whereŷ c (D p θ )) is the prediction from the surrogate model. Then, considering the MLEsβ andφ as fixed values, the posterior meanθ of the calibration parameters θ is deduced according to the prior information and observation data. As pointed out by [8], despite the maximum likelihood estimate plug-in being only approximately Bayesian, the resulting answers seem to be close to those from a full Bayesian analysis.

Sequential D-Optimal Design for Calibration Parameters
How to design the physical and computer experiments efficiently is critical to calibrating the computer models. Usually, the design for the physical and computer experiments is conducted separately in practice. A space-filling design is oftentimes used as the initial design for computer experiments, and some uniform designs or factorial designs are employed for physical experiments. Due to the limitation of physical experiments, after the initial designs, only computer experiments are considered to improve the estimation of calibration parameters sequentially. Here, we propose a sequential D-optimal design for computer experiments, which improves the estimation by maximizing the determinant of the FIM.

D-Optimal Criterion
Following the model presented in Section 2, full output vector d is normally distributed given (θ, β, ϕ). Thus, the FIM can be derived as where p(d|θ, β, ϕ) is the conditional distribution density function, also known as the likelihood function. Similar to [15,21], the formula of FIM of the calibration parameter is presented in the following lemma. The process of the derivation is shown in Appendix A. β, ϕ)). Then, the (i, j)th element of the FIM is where tr(·) is the trace of the matrix, and i, j ∈ {1, · · · , h}.
The corresponding D-optimal design is generated by maximizing which is similar to [22,23], where X denotes the experimental space, and | · | is the determinant of a matrix. The inverse of FIM is the asymptotic variance matrix of the estimate of calibration parameters. Thus, maximizing (10) is equivalent to minimizing the volume of the confidence ellipsoid of the estimate, which can help improve the estimation. Obviously, it is not an easy task to create a D-optimal design by maximizing (10) directly. Here, we utilize the one-point-at-a-time strategy, i.e., add the computer design points by using the Doptimal criterion sequentially, which is referred to as sequential D-optimal design hereafter in this paper.

Algorithm for Generating Sequential D-Optimal Design
Let D p and D c be the initial designs for physical and computer experiments, respectively. In the process of the ith iteration, assume D s = {(x s 1 ,θ 0 ), (x s 2 ,θ 1 ), · · · , (x s i ,θ i−1 )} be the computer experimental points generated sequentially, and y s = {y s 1 , y s 2 , . . . , y s i } be the corresponding computer outputs. Then the estimateθ i can be obtained based on the current full data {(D c , y c ), (D s , y s ), (D p , y p )}. For a fair comparison, the stopping rule is set as 'i ≥ N', where N is the prefixed number of the sequential design points. Then, the next computer experiment point is selected by maximizing Evaluate the computer model at (x s i+1 ,θ i ) and denote the output as y s i+1 . The sequential design is augmented as D s = D s ∪ {(x s i+1 ,θ i )}, and the corresponding computer outputs are augmented as y s = y s ∪ {y s i+1 }. The valueθ N is used as the final estimate of calibration parameter θ. Note that it is an r-dimensional optimization problem to find out the next computer experiment point by maximizing (11), which is a challenge. In order to overcome this problem, discrete optimization is commonly performed. For example, [16] used an algorithm (e.g., [24]) based on a fine grid of the r-dimensional input space. Ref. [17] used a greedy fashion scheme, in which the calibration parameter estimate was considered as the value that maximized the EI criterion over a grid. However, a greedy search based on fine grids may be time-consuming. We employ the following procedure to search for the next design point. Let C be a space-filling design with k points over computer experimental space X and fθ i (c l ) = log |I(θ)| θ=θ i ,x=c l be the corresponding logarithm value of the determinant of FIM at c l andθ i . It is generally reasonable that k = 500 for r < 5 and k = 10 r or 20 r for r ≥ 5, which is similar to [25]. In the numerical simulation studies in Section 4 and the real data analysis in Section 5, since r ≤ 3, then k = 500 is a reasonable choice. Based on C and the corresponding evaluations, we can calculate fθ i (·). Then, the next computer experiment point is selected by maximizing The details about the algorithm to create the sequential D-optimal design are presented as Algorithm 1.

Algorithm 1: Algorithm for Generating Sequential D-optimal Design
Input: Given physical observation data {D p = (x p 1 , . . . , x p q ), y p }, initial computer experiment data {D c = ((x c 1 , t c 1 ), . . . , (x c n , t c n )), y c }. Output: N-run sequential D-optimal design 1 Let D s = ∅ and y s = ∅. 2 Generate a space-filling design C with k points over computer experimental space X . 3 while 0 ≤ i < N do 4 Compute the estimateθ i based on the current full data {(D c , y c ), (D s , y s ), (D p , y p )}. 5 for c l ∈ C do 6 Evaluate fθ i (c l ) = log |I(θ)| θ=θ i ,x=c l .
7 Select x s i+1 by maximizing fθ i (x)|x ∈ C. 8 Evaluate the computer model at (x s i+1 ,θ i ) and denote as y s i+1 .

9
Augment the sequential D-optimal design as D s = D s ∪ {(x s i+1 ,θ i )} and the corresponding outputs as y s = y s ∪ {y s i+1 }.

Simulation Studies
In this section, we investigate the performance of the new proposed sequential Doptimal design using simulation studies. Two numerical simulation examples and one real data analysis are used to compare the performance of the proposed design with the EI [17] and IMSPE designs [16]. To make the comparison fair, the IMSPE design is implemented only regarding computer design in the simulation studies. To evaluate the performance of the designs, the following three statistical metrics are considered: The MSE is used to demonstrate the effectiveness of the estimate of the calibration parameters, which is defined as where ||.|| denotes the corresponding Euclidean distance. Ref. [26] proved that under certain conditions, the Kennedy-O'Hagan calibration estimator converges to the minimizer of the norm of the residual functionδ(x) in the reproducing kernel Hilbert space. As a result, we assume that θ * = arg min ||δ(x)|| N is the best value of calibration parameter. For more details about the reproducing kernel Hilbert space N , please refer to [26,27].θ lm is the estimate of the calibration parameters at the mth replication after l sequential design points, and M is the number of the simulation replications. The MPD is considered to assess the predictive performance of the calibrated computer models, which defined as ||ŷ c (D p ,θ lm ) − y p (D p )|| 2 , for l = 1, · · · , N, ) is the prediction of the computer model withθ lm being the estimate of the calibration parameters at the mth replication with l sequential design points. The MSPE is used to assess the accuracy of the predictions by combining the calibrated computer model and discrepancy, which is defined as whereŷ p (D p ,θ lm ) =ŷ c (D p ,θ lm ) +δ(D p ) is the prediction associated with the physical experiments.

Case Study I
In this section, an example with one calibration parameter and one control variable is considered, which is formulated as 3]. The best value of the calibration parameter is θ * = 1.5 in this case. A constant mean β and a product-form Matérn correlation function with ν = 1 2 are selected as the prior for y c (x, θ). For the calibration parameter, we select the prior of θ to be θ ∼ N(2, 0.01). The size of the physical experiment design (uniform design) is set as q = 10. An MmLHD in [−5, 5] × [0, 3] with 10 points is generated for the initial computer experiment design, i.e., n = 10, and the N = 30 points are to be added sequentially according to Algorithm 1. A total of 100 simulations are performed to calculate the metrics of performance. The results of MSE l are shown in Figure 1. As the increase in computer experimental points sequentially, the calibration parameter approaches the best value, and the proposed method outperforms the other two designs. The results of MPD l are shown in Figure 2, which shares a similar trend with MSE l . The comparison between the original computer model and the calibrated computer model is shown in Figure 3. We can see that the calibrated computer models are closer to the physical observations. As the number of sequential computer experiments increases, the differences between the computer model and physical observations decrease, and the proposed method performs better than the other two designs.   The results of MSPE l are summarized in Table 1, which also shares a similar trend with MPD l and MSE l . The performance of prediction by combining calibrated computer models and discrepancy function is shown in Figure 4. In Figure 4, the predictions by using our method approximate the physical observations best.

Case Study II
In this section, another example which includes three calibration parameters and two control variables and is given in [28] is considered. In this case study, a discrepancy between computer model outputs and physical outputs is also regarding. The physical process is described as The best value of the calibration parameter is θ * = (0.2, 0.3, 0.8). Assume that the expectation of y c (x, θ) is β, and the correlation function of y c (x, θ) is Gaussian. The expectation of δ(x) is set to be zero, and the correlation function of δ(x) is also assumed to be Gaussian. Let the prior for θ 1 , θ 2 , and θ 3 be U(0, 0.25), U(0, 0.5), and U(0, 1), respectively. A total of 500 simulations are performed to assess the performances of the sequential designs. A design with 10 points are employed for physical experiments, and a MmLHD with n = 20 points over [−1, 1] 3 × [0, 0.25] × [0, 0.5] is utilized as the initial design for computer experiments. In this case, the number of sequential points is set to be N = 40. The results of MSE l are drawn in Figure 5. As the size of sequential experiment points increases, the calibration parameter approaches the best value. The proposed method has smaller MSEs than the other two methods, which means the proposed method performs better. Figure 6 shows the results of MPD l , which also demonstrate the effectiveness of the proposed design for calibration.  The simulation results in Figures 1 and 5 show that, regardless of whether the discrepancy function exists or not, with the increase in sequential points, the estimate of calibration parameters gradually approaches its best value, indicating that the estimation of calibration parameters converges. However, the strict mathematical proof of this conclusion is a complicated problem, which cannot be solved in this paper and will be paid attention to in a future study.
From Figure 6, we can easily find out that when 20 points are sequentially added, the discrepancy tends to be stable. As shown in Figure 7, the physical prediction combining discrepancy and calibrated computer model by utilizing the proposed sequential D-optimal method is the closest to the physical observations, followed by the IMSPE method and EI method. The results regarding the MSPE are presented in Table 2, which shares similar conclusions as those in the case study I.

Real Data Analysis
In this section, a real data example with three control inputs and one calibration input is considered, which is presented in [8]. Figure 8 shows a concise description of the resistance spot welding process. Two metal sheets of a particular thickness (thickness) are compressed through two electrodes under a specifically applied load (load). A direct current of a certain magnitude (current) passes through the sheets via the two electrodes, and the heat produced by the current flow causes the welding surfaces to melt. After cooling, a weld nugget with a specific dimension (diameter) is formed, which is of particular interest. In this manner, the two metal plates are welded. The resistance at the contact surface is particularly critical in determining the magnitude of heat generated. Because the contact resistance at the contact surface is not well understood as a function of temperature, the calibration parameter is specified and adjusted based on the field data. The effect of this calibration parameter on the behavior of the model is a focus in this case. Ref. [8] comprehensively described the inputs for this example. According to the evaluation of the model developer, the verification experiment focuses on three control inputs (thickness, load, and current). Table 3 lists the control and calibration inputs and the corresponding intervals. A constant mean β and the product-form exponential correlation function are considered as the prior for y c (x, θ). As (2) shows, a zero mean and the Gaussian correlation function are considered as the prior for δ(x), and more details can be found in [2,6,28,29], etc. We randomly choose 10 physical experiments from the non-replicated 12 physical experiments. An initial computer design with 20 points is generated and the corresponding outputs are simulated according to [8]. N = 40 points are added sequentially according to Algorithm 1. The results about the MSPE are summarized in Table 4, which illustrates the superior of the new proposed method due to the smaller MSPE values.

Conclusions and Remark
In order to reduce the cost of physical experiments, mathematical models or computer models are utilized to approximate the real physical process. However, the computer models' fidelity to physical process depends on the physical unobservable calibration parameters. This paper proposed a D-optimal design to augment computer experiments sequentially by regarding the Kennedy and O'Hagan model. By using the D-optimal criterion, computer design points are selected sequentially to gather more comprehensive information to reduce the uncertainty about the estimate of the calibration parameter. Then, the computer model can mimic the real physical process well with the tuned calibration parameter. Simulation studies are made to assess the performance of the newly proposed method compared with EI and IMSPE methods. The results show that the newly proposed method outperforms the other two methods in terms of MSE, MPD, and MSPE. An analysis based on the real data introduced in [8] also demonstrates the superior performance of the new method.