Iterative Method for Tuning Complex Simulation Code

Tuning a complex simulation code refers to the process of improving the agreement of a code calculation with respect to a set of experimental data by adjusting parameters implemented in the code. This process belongs to the class of inverse problems or model calibration. For this problem, the approximated nonlinear least squares (ANLS) method based on a Gaussian process (GP) metamodel has been employed by some researchers. A potential drawback of the ANLS method is that the metamodel is built only once and not updated thereafter. To address this difficulty, we propose an iterative algorithm in this study. In the proposed algorithm, the parameters of the simulation code and GP metamodel are alternatively re-estimated and updated by maximum likelihood estimation and the ANLS method. This algorithm uses both computer and experimental data repeatedly until convergence. A study using toy-models including inexact computer code with bias terms reveals that the proposed algorithm performs better than the ANLS method and the conditional-likelihood-based approach. Finally, an application to a nuclear fusion simulation code is illustrated.


Introduction
Modern computer simulation codes contain various unknown parameters. Assuming the validity of the simulation code, we can adjust or estimate such parameters using the nonlinear least squares estimation (NLSE) method, which minimizes the sum of the squared differences between computer responses and real observations. This procedure is called calibration (e.g., Kennedy and O'Hagan 2001; Higdon et al. 2008; Tuo and Wu 2018) or code tuning (e.g., Cox, Park, and Singer 2001;Kumar 2015). It is formally defined as the process of improving the agreement of a code calculation or set of code calculations with respect to a chosen and fixed set of experimental data via adjustment of the parameters implemented in the code (Trucano et al. 2006). Han, Santner, and Rawlinson (2009) differentiated between tuning parameter and calibration parameter. In this study, however, the two parameters are treated as the same, henceforth it is referred to as the tuning parameter.
If a simulation program is complex, with one execution requiring several hours, the NLSE method may not be computationally feasible. In this case, a statistical metamodel can be built to approximate the unknown functional relationship between a set of controllable input variables and a simulated response variable. This metamodel is then employed in place of the original simulation code in the NLSE method, making the problem solvable and computationally feasible. This method was described by Cox, Park, and Singer (2001), where a Gaussian process (GP) model was employed as the metamodel of a complex simulation code. It is called the approximated NLS (ANLS) method. As alternatives to the ANLS, Cox, Park, and Singer (2001) proposed likelihood-based methods. Kennedy  A potential drawback of the ANLS method is that the metamodel is built only once and not updated thereafter; thus, the computer data are no longer used live. To address this, in the present report, an iterative algorithm is proposed, in which the likelihood function is maximized, and the squared distance is minimized iteratively until convergence is achieved. That is, the tuning parameters of the simulation code and model parameters of the GP are repeatedly reestimated and updated. The parameters of the GP model are estimated from updated combined data using the maximum likelihood approach. Then, the tuning parameters are re-estimated using the ANLS method. These two optimizations are iteratively executed until convergence is achieved. We call this method the "Max-min" algorithm.
The remainder of this report proceeds as follows. Section 2 presents a GP model to approximate a complex simulation code and real experimental data. Section 3 describes the calibration methods, including our proposed method, and the calculation of standard errors. Section 4 presents a toy-model simulation study. Section 5 applies the proposed methods to a nuclear fusion simulator, and Section 6 summarizes the results and provides suggestions for future research. Some details of this study are provided in the Supplemental Material.

Gaussian process model as a surrogate
For the metamodel of the simulation code, we use a GP or a spatial regression model that treats response y(x) as a realization of a random function superimposed on a regression model: where f s are known functions, and βs are unknown regression coefficients. Here, the random process Z(.), which represents a departure from the assumed linear model, is assumed to be a GP with mean zero and covariance cov(t, u) = V (t, u) = σ 2 R(t, u) between Z(t) and Z(u) for t = (t 1 , ..., t d ), u = (u 1 , ..., u d ), where σ 2 represents the process variance (a scale factor), and R(t, u) is the correlation function. When the response of a computer code is stochastic, the random component term is added to the model (1). However, we do not include in this study because a computer code is assumed to be deterministic. Some possible choices for the correlation function are obtained from the Gaussian correlations denoted by where θ ≥ 0. These are special cases of a power exponential family with a power of 2. The non-negative parameter θ determines the covariance structure of Z: a small θ reflects high correlations between nearby observations, whereas a large θ reflects low nearby correlations. One may consider a different version of (3) by taking several θ values as follows: where θ i ≥ 0, for i = 1, 2, · · · , d. This is a legitimate separable correlation function because it is a product of valid correlation functions. We refer the readers to Santner, Williams, and Notz (2018) for more information on this GP model and its application to the design and analysis of computer experiments. Once the data have been collected at the design sites or at "training" inputs X, the parameters are estimated via the maximum likelihood estimation (MLE) method and are then plugged in to predict y(x 0 ) as in (5), where x 0 is an untried site or a "test" input. This prediction is called Kriging. The empirical best linear unbiased prediction  with the MLEs of the parameters is denoted bŷ where f 0 is the known linear regression function vector, F is a design matrix, r t 0 is the correlation vector between Y (x 0 ) and model outputs Y (X), y is the vector of observations collected at the design sites, andβ is the generalized least squares estimator of β (see the Supplemental Material for the details).
The combinations of βs and θs determine the model, but the following simple GP model is considered first: with the correlation function (3) or (4). When a common θ of (3) is used, we call (6) "Model 1" in this study. When the several θ i s of (4) are used, we call (6) "Model 2".
3 Methods for code tuning

Data structure
For notational convenience, experimental data is denoted by the subscript "E," computer simulation data by the subscript "C," and "both" computer and experimental data by "B." Let τ be an adjustable parameter vector to be estimated. Let T be the input variables of the computer code corresponding to τ . The original experimental input variables are denoted by X. Let q and p be the dimensions of τ and X. Further, let n C andn E be the number of observations; then, n B = n C + n E . The details of the data structure for calibration are described in the Supplemental Material.

Approximate nonlinear least squares
In this subsection, the GP model approach to be used in the ANLS procedure provided in Cox, Park, and Singer (2001) is described. If the parameters σ 2 , β, and θ are known, then, for a given value of τ , a "prediction" of y E (τ, x iE ) can be calculated for the given computer data. This is obtained using equation (5) and X E , F E , y E , and V EE . Here, the computer data alone are used to calculate β andθ, and the data are not used thereafter. Note that X E , F E , R(X E , X C ), R(X E , X E ), and β are functions of τ . A design site selected for a computer experiment is denoted by (T , x C ). Then, the computer response y (or y C ) at (T , x C ) is where Y represents the expected value of the output from computer code. Y can be viewed as the value produced from a theory. Since we assume that the computer code is close to the real experimental data with some variation if the tuning parameters in the computer code are optimal, the response of the real experiment y E at (τ , x E ) is modeled by Here Y (τ , x E ) also represents the expected value of the response in the real experiment. This common Y in the two abovementioned equations (7) and (8) connects the computer code and the real experiment. The stochastic term E is assumed to be independent and identically distributed with mean zero and variance σ 2 E . On the other hand, a computer code is sometimes considered inexact because of computer model bias, which is the systematic difference between the model and the truth (Kennedy and O'Hagan 2001; Gramacy 2016; Plumlee 2017). In that case, the response of the real experiment can be modeled by where b(τ , x E ) represents the model bias or discrepancy. A drawback of the computer experiment is that one run of a complex simulator sometimes requires several minutes, say, m minutes. Then, τ is usually estimated by minimizing the residual sum of the squares: where y E i is an observed response from the real experiments, and Y (τ , x E i ) is the expected value of the output from the computer code at the experimental point (τ , x E i ). One evaluation of RSS(τ ) requires approximately n E × m minutes. It is thus computationally infeasible to run the code as many times as needed for an iterative nonlinear optimizer to find τ . The ANLS method first fits the GP Model 1 or Model 2 as defined in equation (6), with MLE using computer data alone. Then, upon treating the fitted prediction model as if it were the true model, the method attempts to determine the τ that minimizes the residual sum of squares with predictors: whereŶ (τ , x E i ) is the empirical best linear unbiased prediction of Y (τ , x E i ), as in (5). Since it is difficult to have a closed-form minimizer of (11), a numerical optimization routine is necessary to determine τ . Note thatŶ is a computationally cheaper emulator (surrogate or metamodel) of the expensive simulation code. This makes the problem computationally feasible. It is applied for each GP model (Model 1 and Model 2) stated in (6). The advantages of this method are that it is reasonably economical and easy to implement, and the computer and experimental data are uncoupled. The ANLS method does not require the functional relation between the inputs and output to be the same for the computer code and real experiments. In contrast, the likelihood-based approaches described in the next subsection do require this, but they employ marginal likelihoods to estimate model parameters. The prediction residuals, r i = y E i −Ŷ ( τ , x E i ), can be used to check the validity of prediction model, τ , and of the ANLS method. A potential drawback of the ANLS method is that it does not account for uncertainty in the approximation of Y C byŶ . Another difficulty is that the metamodelŶ is built only once and is not updated thereafter.

Likelihood-based tuning methods
Given a computer code and our GP approach for y, a unified statistical approach is available. We have the likelihood for all the parameters, including the tuning parameters τ ; the error term parameter σ 2 E , and the random function parameters β, θ, and σ 2 . Thus, all parameters can be estimated using the MLE method. We refer to this method as the full MLE. It can be done for each GP model defined in (6). The −2 times concentrated log likelihood function (except for constants) of all parameters for the combined data withβ B andσ 2 B plugged in is where Some other approaches based on the likelihood function, including the full MLE, were also proposed by Cox, Park, and Singer (2001). One of them is the "Separated MLE" (SMLE) method, which maximizes the conditional likelihood function of experimental data when the computer data is given. Here, the parameters β, θ, and σ 2 are estimated by maximizing the marginal likelihood for the computer data only. These are plugged into the conditional likelihood of the experimental data given the computer data. This is then maximized with respect to γ E and τ to obtain estimates of those parameters.
One advantage of the likelihood-based method is that it can simultaneously use both computer and experimental data to estimate τ , whereas ANLS method uses computer data only. These likelihood-based approaches enrich the tuning methods. Cox, Park, and Singer (2001), found the SMLE method to be better than the full MLE method. Thus, SMLE is compared with the proposed method in this study. The details of the SMLE are provided in the Supplemental Material.

Proposed method 4.1 Max-min algorithm
The following are the steps for the proposed tuning method. We call it a Max-min algorithm because it uses maximization and minimization iteratively.

Algorithm Max-min: iterate Step 3 and Step 4 until convergence is achieved
Step 1 (model building): build a surrogate (6) using the MLE for the given computer data only.
Step 3 (maximization): build a new surrogate (6), using the MLE for the combined data with the fixed τ obtained in the previous step.
Step 4 (minimization): set iteration i = i + 1, and find τ by minimizing RSS p (τ ) in (11) using the surrogate built in Step 3. If τ satisfies the stopping rule, then stop; otherwise, go to Step 3. Note that in each iteration of Steps 3 and 4, τ is updated; thus, the estimates of the parameters of θ, β, σ 2 E , and σ 2 are updated. We expect this to positively influence the finding of τ of Step 4. Steps 2 and 4 are the same in terms of minimizing RSS p (τ ), but Step 2 uses only computer data, while Step 4 uses either the combined or computer data. Steps 1 and 3 are the same in terms of obtaining the MLE by maximizing the likelihood function, but Step 1 uses only computer data, while Step 3 uses the combined data. The likelihood function in Step 3 is L(η B ; τ , y B , X B ), which is used for the full MLE method. But here, the tuning parameters are fixed as the τ that was obtained in the previous step. For the optimizations in Steps 2, 3, and 4, quasi-Newton numerical algorithms were used.
On using the combined data in Step 3, we assume that the functional relation between the inputs and output is the same for the computer code and physical process. The computer data and the experimental data are linked by the use of this common response function, with τ being the value for the tuning parameters in data from the physical process.
The Max-min algorithm stops when one of the following rules is satisfied: for i = 1, 2, · · · , 1 (maximum iteration): Number of iterations reaches the pre-assigned maximum number, where 'ftol' is a pre-assigned small value for tolerance. Here RSS p ( τ 1 ) is the minimum value of RSS p obtained in Step 2, and RSS p ( τ i ) is the minimum value obtained in Step 4 in the i-th iteration.
When the RSS p in Step 4 is greater than that of Step 2 or that of the last iteration, a small random fluctuation on τ is given before Step 3. Without this fluctuation setting, the algorithm stopped within four iterations. Based on our experience, the random fluctuation caused a reduction in RSS p . However, it did not make the algorithm execute more than 20 iterations.
When x 0 is a given site representing ( τ , x E i ) in the experimental data (X E ), the prediction formula needed for computing RSS p in Step 2 iŝ where f 0 is the known linear regression function vector; r 0C is the n C × 1 correlation vector between Y (x 0 ) and Y (X C ), and The prediction formula needed for computing RSS p in Step 4 iŝ where r 0B is the (n E + n C ) × 1 correlation vector between Y (x 0 ) and Y (X B ). Here,θ B , β B , andγ E are the MLEs from the combined data in Step 3. Note thatθ B andγ E are needed for constructing V B . Even though x 0 and X E are included in the construction of r 0B and V B , the predictionŶ (x 0 ) in (16) is not an exact interpolation because of the positiveγ E . Our approach was motivated by the iteratively re-weighted least squares method in regression analysis. The alternating estimation of parameters recursively is similar to the EM (expectationmaximization) algorithm. One can view this method as similar to a frequentist version of the Bayesian modularization approach (Liu, Bayarri, and Berger 2009). Some practical suggestions for the calibration of large-scale simulations and optimization are provided in Gramacy (2016, Section 4). However, our study was conducted independently from their study.
The advantages of the proposed Max-min algorithm are as follows. The uncertainty in the approximation of Y usingŶ in Step 4 is smaller than that in the ordinary ANLS, because the MLE method in Step 3 and the prediction model in Step 4 use combined data (a larger sample size), while the ordinary ANLS method uses only the computer data in Step 1 for MLE and in Step 2 for prediction. In addition, the Max-min algorithm accounts for the covariance between the computer data and experimental data in building the metamodel. Moreover, our test function experience revealed that the solutions from the Max-min algorithm are less influenced by the initial value of τ than by the ANLS method because the estimated τ in the Max-min is iteratively updated several times.

Approximated confidence region of estimates
Once τ has been estimated, some indications of the accuracy of the estimates are generally necessary. Thus, we rely on the asymptotic theory for the nonlinear least squares estimator that appears in the regression analysis (Draper and Smith 1981). The approximate 100(1 − α)% confidence region of τ is obtained by where F α (q, n E −q) is the upper (1−α) percentile of the F distribution with q and n E −q degrees of freedom (q is the number of parameters in τ ). Wong, Storlie and Lee (2017) considered a Bootstrap approach to compute the uncertainty of the estimates. 5 Toy-model study 5

.1 Exact computer models
In this section, we apply the methodology outlined in the previous sections to the seven test functions Y (τ , x) that are easy to compute. In the first five test functions, the experimental data with sample size n E are generated by where the random variable e follows a normal distribution with mean zero and variance σ 2 e . These five examples assume that the functional relation between the output and inputs is the same for both the computer code and the physical process. By contrast, in the last two test functions (6 and 7), the model bias term b(x) is added to the function Y (τ * , x) in (18). Here, τ * stands for the true tuning parameters and is specified for each toy model. No random error is given to the computer responses. Sample sizes for both data are set to 30 (n C = n E = 30) for the first five test functions. For generating the computer and experimental data, the same uniform distributions were used for the random numbers of x variables. The details on the eight test functions are described in the Supplemental Material.
To address the uncertainty resulting from the designs for the input variables, we repeated the τ estimation using 30 random Latin-hypercube designs. Thus, 30 sets of estimates of τ were obtained, and the averages and standard deviations are reported. Two different GP models of (6) were used with the correlation functions of (3) and (4), which are called "Model 1" and "Model 2," respectively. Figure 1 shows box plots of the distance to the true value in five test functions in which Model 2 is employed as a surrogate. The Max-min algorithm generally works better than the ANLS and SMLE methods. Tables S1 to S8 in the Supplemental Material present the results for each toy model, each method, and Model 1 and Model 2. Some columns show the average of Euclidian distances (Dist) between the true τ and estimates, with standard deviations in parentheses. The last column shows the mean squared error (MSE) of the estimates obtained using the following formula: where std(τ i ) is the standard deviation of each estimate computed from 30 repetitions. Figures S1 to S5 show the box plots for each test function and method. In terms of the average distance to the true value and MSE in Tables S1 to S5, the Max-min algorithm works better than ANLS for all toy models. Max-min also works better than SMLE for test functions 1, 2, and 4. Finally, SMLE works slightly better than ANLS, except in test function 2.    Figure S6 shows the typical convergence of the Max-min algorithm for 10 trials of test function 1. In most cases, the algorithm stops at the third or fourth iteration, except for a few cases wherein it stops at the sixth or seventh iteration. This means that the improvement of τ in the second iteration is significant, while it may not be significant after the second iteration. The small random fluctuation on τ , as mentioned in subsection 4.1, was not applied in this computation.

Another prediction
A potential drawback of the Max-min algorithm is that it partially fails to decouple between the estimation of τ and building of the metamodel because it uses the combined data in both Steps 3 and 4. To address this difficulty, instead of using equation (16) in Step 4, one can employ prediction (15) but withθ B and β B (or β C ): This leads to another version of the Max-min algorithm, which is applied to the next toy models. Note that this version does not useγ E . We denote this prediction asŶ C|B (x 0 ). We can employ β C instead of β B in (16), but we have not yet attempted to do so.

Inexact computer models
A computer model is often considered inexact. This means that the computer model does not perfectly match the real system even if some parameters included with it are optimal (Plumlee 2017). This is the computer model bias; it is the difference between the model and the truth. The model-bias term b(x) is included in the following two test functions. These are employed to check the performance of the proposed method for cases where the functions generating the output for the computer code and for the physical experiment differ, but we incorrectly assume that (7) and (8) hold.
Test function 6: n C = n E = 20 Test function 7: n C = n E = 20 Computer data : Test function 6 is modified from Plumlee (2017). Test function 7 was used in Bastos and O'Hagan (2009), in Goh et al.(2013), and in Gramacy (2016). Twenty random Latin hypercube designs were used repeatedly to address the uncertainty.
In the previous exact toy-model study, we used the difference between τ and true τ for the physical process as a measure of performance. This requires that there be a true τ , but there would be no such true τ in an inexact computer model. In this situation, the minimum value of RSS p would be a more appropriate measure.
In calculating RSS p in the Max-min algorithm for the above two test funtions, we set ftol = 1.e-4, and maxagain = 7. For the small fluctuation of τ , random numbers from N (0, σ 2 τ ) were used where σ τ = max( τ × 0.1, 0.3). Figure 2 shows parallel coordinated box plots of RSS p values for test function 6 computed from 20 Latin hypercube designs. RSS p values are calculated using the ANLS method and the Max-min algorithm with GP Model 1 and Model 2 via two different predictions. The predictions byŶ B in (16) andŶ C|B in (20) in calculating RSS p in the Max-min algorithm were attempted for the sake of comparison. The Max-min algorithm worked better than ANLS in most cases. The left panel of Figure 3 is an ANOVA-type main-effect plot based on mean values, which shows improvements from ANLS to Max-min, from Model 1 to Model 2, and fromŶ C|B toŶ B . The mean values from the last two predictions (IB0 and IB1) are calculated within the Max-min algorithm.
Relative improvement (RI) from ANLS to Max-min is provided in Table 1. RI is defined by     For tuning test function 7, we considered a simple bias correction (BC) technique. Two constants for additive and multiplicative corrections were set for the predictor (Fernandez-Godino et al. 2016) asŶ This naive BC method may not be better than sophisticated methods such as that using ρ(x) and δ(x), or the approach by Kennedy and O'Hagan (2001), but we hope for this to work better   Figure 4: This is the same as Figure 2, but for test function 7, calculated without (non-BC) and with bias correction (BC).
than the non-BC case. The objective function RSS p in (11) is now changed to Then, τ , ρ, and δ are estimated simultaneously by minimizing RSS bc p in both the ANLS and the Max-min methods. Figure 4 shows parallel coordinated box plots of the RSS p values for test function 7, computed from 20 Latin hypercube designs. The left panel is for non-BC, while the right one is for BC. The acronyms in this figure are the same as in Figure 2. This figure shows that the Max-min algorithm worked better than ANLS in most cases. The right panel of Figure 3 is a main-effect plot for test function 7. It shows improvements from ANLS to Max-min, Model 1 to Model 2, non-BC to BC, andŶ C|B toŶ B . Note that the patterns of the main effects for test functions 6 and 7 are very similar. The RI from ANLS to Max-min is provided in Table 1. The bias correction was most effective when used with Model 2 and with the predictionŶ B .
6 Application to nuclear fusion model 6

.1 Nuclear fusion data
A simple measure of energy efficiency in a nuclear fusion device (called a tokamak, from the Russian language) is the global energy confinement time τ E . The theoretically based confinement model can be written as follows (Kay and Goldston 1985): P, I, N, B) , (24) where f is a known function calculated using a complex simulation code called Baldur, P is the total input power, I is the plasma current, N is the electron density, B is the magnetic field, and τ = ( τ 1 , τ 2 , τ 3 , τ 4 ) are the following adjustable parameters that determine energy transfer, that is, drift waves, rippling, resistive ballooning, and the critical value of η i (which provokes increased ion energy losses for the drift waves), respectively. The experimental data comprises only P, I, N, B, and the real observation y E , whereas the computer data comprises eight independent variables (T , P, I, N, B) and computer response y C obtained using Baldur. The experimental data were drawn from the database of S. Kaye: 42 observations from the PDX (Poloidal Divertor Experiment) tokamak in Princeton and 64 from the Baldur simulator (Singer et al. 1988). Because the Baldur simulator requires five CPU minutes on a Cray supercomputer for one execution, a careful selection of input points is required, which is a statistical design problem for a complex simulation code. For this purpose, we used a data-adaptive sequential optimal design strategy, which is described briefly in the Supplemental Material.  Table S9 in the Supplemental Material presents the MLE of the parameters of GP Model 1 obtained from the combined data. Table 2 provides the results of τ estimation obtained using ANLS, SMLE, and the Max-min methods on the basis of Model 1. A quasi-Newton optimization routine was employed in searching for τ . Several starting values were tried to avoid the local minima. The last column in Table 2 shows the value of RSS p at the convergence of the algorithms. The RSS p for SMLE was obtained by calculating RSS p for the estimated τ . The RSS p for the Max-min with Model 1 is the smallest among the three methods. Figures 5 and 6 show the residual plot (residual vs. predicted values) and confidence regions of the tuning parameters (τ 1 and τ 4 ) for the tokamak data, which were obtained using the Maxmin algorithm with Model 1. In Figure 5, the predicted values for the computer data (circles) had a wider horizontal range than those for the experimental predicted values (crosses). In addition, the residuals for the computer data had a narrower vertical range than those for the experimental data. These results indicate that the computer data are fitted relatively well, and that good coverage of the range of experimental observations is ensured.

Summary and Discussion
Using a GP model, we considered an iteratively re-estimated ANLS method to tune complex computer code to data, namely a Max-min algorithm. This method is an extension of the ANLS method. A simulation study using toy functions suggests that the proposed Max-min algorithm works better than the ANLS and SMLE methods. We applied this technique to a computational nuclear fusion model. The proposed method can be useful for other applications in disciplines in which unknown theoretical parameters must be estimated using complex computer codes and real experimental data. Moreover, the Max-min algorithm is also applicable to other metamodels, such as spline (Wong, Storlie and Lee 2017), support vector machines, and neural networks. The Max-min algorithm requires more computing time than the ANLS method. This is because Max-min uses the combined data, resulting in a larger correlation matrix that must be inverted. It also requires more iterative minimizations of RSS p compared to the ANLS method, which requires just one iteration. Two versions of predictions were tried in the last two toy models wherein bias functions were added. The prediction based on both data sources (Ŷ B ) reduced RSS p (residual sum of squares of prediction) more than that based on computer data given the parameter estimates obtained from both data sources (Ŷ C|B ). Even thoughŶ B reduced RSS p more thanŶ C|B , the improvement of the estimation of true τ compared toŶ C|B is not known. A simple bias correction was considered in test function 7. This bias correction was the most effective when used with GP Model 2 and with the prediction based on both data sources.
We tried two different GP models in the toy-model study, and found that the estimates of tuning parameters were changed according to the selected GP model. Model 1 in (6) is a superimposition of a GP with correlation function (3) on the first-order regression model. We believe that the first-order regression part sometimes works well to fit the dominant relationship, even though the isotropic correlation (3) may be unrealistic for capturing the variability of responses from the dominant relationship. It is computationally easier to obtain the MLE when the number of input variables is large. Model 2, with (4) in (6) assumes a separable Gaussian correlation, which is a popular family of correlation models as per the literature (Santner, Williams and Notz 2018), but it is computationally expensive to obtain the MLE if many inputs exist.
A model selection procedure among many GP models with various combinations of nonzero β j s and θ j s in (6), as in Marrel et al. (2008), Chen et al. (2016), and Lee and Park (2017), may lead to realization of a better surrogate model. In this study, we did not consider the identifiability of the tuning parameter, which is an important focus of recent developments (e.g., Plumlee 2017; Tuo and Wu 2018).
There are, however, basic limitations with tuning computer code to real-world data regarding the experimental design. We found that the performance of tuning methods is significantly dependent on the designs for both the computer experiments and the physical experiments. Some authors, including Cailliez, Bourasseau, and Pernot (2014), and Beck and Guillas (2016), explored this topic. We agree that a sequential tuning approach is practically useful, as in Pratola et al. (2013), Kumar (2015), and Damblin et al. (2018). Further research on relevant designs under the sequential tuning approach will be helpful.

A.1 Maximum likelihood estimation in Gaussian process model
Once the data have been collecte at the observation sites {x 1 , ..., x n }, we use the maximum likelihood estimation (MLE) method to estimate the parameters in linear model part and covariance function. Since we assume that y(x) is a Gaussian process with mean F β and covariance matrix σ 2 R, the likelihood function of y is where F is a so-called design matrix. When the covariance parameters θ and γ C are specified, the MLEs of σ 2 and β are denoted by Here,β is the generalized least squares estimator of β. Since the likelihood equations do not lead to a closed-form solution, a numerical optimization procedure is required. The Cholesky decomposition V = U t U is used as a major computation in calculating the likelihood function, where U is an upper triangular Cholesky factor. The computational details of calculating and minimizing negative log-likelihood function are provided in Park and Baek (2001). One can use a R program DiceKrig (Roustant 2012).

A.2.1 Computer and experimental data
For notational convenience, experimental data is denoted by the subscript "E" and computer simulation data by subscript "C." Let τ be an adjustable parameter vector to be estimated. Let T be the input variables of the computer code corresponding to τ . Here, τ is a vector of the deterministic tuning parameters and T is a vector of the random variables. The original experimental input variables is denoted by X. Let q and p be the dimensions of τ and X. Further, let n C , n E be the number of observations. Then, we have the data matrix of the independent variables: X C and X E for computer and experimental data; Here, t ij in X C represents the j-th value of the i-th T variable (T i ) and x Eij and x Cij denote the j-th value of the i-th X variable of experimental (X Ei ) and computer (X Ci ) input. X E is a n E × (q + p) matrix and X C is a n C × (q + p) matrix. Note that the first part of X E is composed of the unknown parameters τ 1 , · · · , τ q , while the corresponding part of X C comprises input values (t ij ).

A.2.2 Combined data
The following notations for combined computer and experimental data are introduced: for the data matrix of the independent variables; the so-called "design matrix, defined as the functions of the values of input variables; and the computer responses and real observations, respectively. Here, the subscript "B" indicates the combined "both" computer and experimental data. Note that variables T are incorporated in the simulation code as design sites. X C and F C contain T , while X E and F E are the functions of the unknown parameters τ . The Gaussian process model is now simultaneously applied to computer and experimental data. Let η = (τ , θ, γ C , γ E , σ 2 , β), where γ C = σ 2 C /σ 2 and γ E = σ 2 E /σ 2 , which are the variance ratios for the computer and experimental data. Here, σ 2 C and σ 2 E are the variances of error term ( ) in the Gaussian process model for the computer and experimental data, respectively. When necessary, β C and β E are used to denote the regression coefficients for the computer and real experimental data. Then, given the independence and normality assumptions, we have where where R(X C , X E ) represents a n C × n E matrix composed of the correlations computed between X C and X E . Note that V B is a n B × n B positive definite covariance matrix for the combined data, where n B = n C + n E . We set γ C = 0 because only a deterministic computer model is considered in this study.

A.3 Separated MLE for code tuning
For the details of SMLE, we make use of the conditional distribution of the experimental data given the computer data, which is normally distributed with mean and covariance where covariance matrices V 's are given as in (32). In these formulae, we suppressed the parameter dependencies in µ E = F E β E , µ C = F C β C , V CE , V CC , and V EE . Now the −2 times concentrated log conditional likelihood function (except for constants) withβ andσ 2 E|C plugged in is −2 log L(τ , γ E ; X E , y E | y C ,γ C , β C ,θ, X C ) = n E logσ 2 E|C + log |V E|C |, A. 4 Test functions for toy model study Test function 1: nC = nE = 30 for all five test functions.
In each test function, n C values of (T , x) were selected as inputs for the "computer code", that is, the function Y (T , x) is evaluated at n C values. The inputs were chosen to be well spread around a reasonable space known to potentially contain the true parameter value. n C computer data points and n E experimental data points were generated by using random Latin-hypercube designs, except the τ values were used instead of T . In real application situation, one would consider more sophisticated designs such as data-adaptive sequential optimal experiments as reported in the next section.

A.5 Sequential designs in tuning a nuclear fusion simulator
For given a Gaussian process model, the A-optimal design is obtained by minimizing the integrated mean squared error of prediction (MSEP) with respect to a design S , where Q is the design region, and µ is a "weight function" which may be the empirical measure of uniformly distributed random points. Note that neither M SEP nor IM SE depend on the unknown parameters β and σ 2 , but depend on θ, γ C and design S. This makes it possible to design an experiment (for specified values of θ and γ C ) before taking the data. Because θ is generally not available for the initial design stage, in our example, we used an rough estimate of θ based on a previous similar work given by a Baldur specialist. Initially we found 10 optimal design points for eight variables (T , P, I, N, B) which minimize IM SE over the design region Q, with θ = .5 and γ C = .001.
The following is the data-adaptive sequential optimal design procedure that we have used in this study.
Step 1. Collect computer observations based on the given optimal design (n 1 points, say).
Step 2. Find an appropriate model and estimates of parameters (θ, β, σ 2 , γ C ) from the computer data.
Step 3. Check the MMSE, and stop constructing the next stage design if MMSE is smaller than a preassigned target value. Otherwise, go to the next step.
Step 4. Use the estimates and model found in Step 2 to choose the next stage optimal design (n 2 points, say) under the condition that the previous design is given (i.e, update n 2 more points to the previous design to make n 1 + n 2 points).
Step 5. Collect (n 2 ) more observations, and go to Step 2.
The MMSE (maximum mean squared error of prediction) is defined as where x i , i = 1, 2, · · · , K, are d-dimensional random vectors. Note that the MMSE is used here as a measure of accuracy of a given prediction model.        Regression coefficient for T 1 -0.027 β 2 Regression coefficient for T 2 -0.010 β 3 Regression coefficient for T 3 0.001 β 4 Regression coefficient for T 4 -0.015 β 5 Regression coefficient for P -0.031 β 6 Regression coefficient for I 0.009 β 7 Regression coefficient for N -0.012 β 8 Regression coefficient for B 0.017 σ 2 Variance of Y 1.506E-04 γ E Variance of σ 2 E /σ 2 0.454  Figure S1 but the test function 2.  Figure 9: Same as Figure S1 but the test function 3.  Figure 11: Same as Figure S1 but the test function 5.