Abstract
While most calibration methods focus on inferring a set of model parameters that are unknown but assumed to be constant, many models have parameters that have a functional relation with the controllable input variables. Formulating a low-dimensional approximation of these calibration functions allows modelers to use low-fidelity models to explore phenomena at lengths and time scales unattainable with their high-fidelity sources. While functional calibration methods are available for low-dimensional problems (e.g., one to three unknown calibration functions), exploring high-dimensional spaces of unknown calibration functions (e.g., more than ten) is still a challenging task due to its computational cost and the risk for identifiability issues. To address this challenge, we introduce a semiparametric calibration method that uses an approximate Bayesian computation scheme to quantify the uncertainty in the unknown calibration functions and uses this insight to identify what functions can be replaced with low-dimensional approximations. Through a test problem and a coarse-grained model of an epoxy resin, we demonstrate that the introduced method enables the identification of a low-dimensional set of calibration functions with a limited compromise in calibration accuracy. The novelty of the presented method is the ability to synthesize domain knowledge from various sources (i.e., physical experiments, simulation models, and expert insight) to enable high-dimensional functional calibration without the need for prior knowledge on the class of unknown calibration functions.
1 Introduction
Simulation models provide a low-fidelity approximation of physical systems but need to be thoroughly validated before being used to inform design decisions. For example, a high-fidelity response model as given by the red curve in Fig. 1 can be approximated through a low-fidelity response model as given by the blue surface. However, low-fidelity models typically include an additional set of model parameters t that need to be calibrated using a set of low-fidelity z and high-fidelity y data points [1–5]. These additional parameters are often introduced using approximation models (e.g., force field parameters associated with coarse-grained beads in molecular dynamics) or phenomenological models (e.g., constitutive law parameters in a finite element simulation). However, in many engineering challenges the model parameterst have a functional relation with the model variablesx (i.e., as represented by the yellow curve on the horizontal plane). For example, such challenges are encountered when calibrating the degradation of industrial equipment health indicators as a function of time [6,7], the calibration of interphase properties in composite materials as a function of the local structure [8], and the calibration of force field parameters in a coarse-grained polymer model with respect to temperature [9]. While a piece-wise linear approximation of the calibration functions can be achieved [9], this would also introduce a large number of parameters and typically results in unreasonably complex functions. Specifically, these piece-wise functions are likely to lead to erroneous results when evaluating the low-fidelity model with different initial conditions or length scales (e.g., a calibrated coarse-grained (CG) molecular dynamics simulation can be used to study physical phenomena at length and timescales that are too computationally exhaustive to be studied using the high-fidelity all-atomistic (AA) model). Consequently, we would therefore desire a statistical method to systematically balance the low-fidelity model complexity (i.e., reduce the number of model parameters) with the calibration accuracy (i.e., the low-fidelity model and high-fidelity model discrepancy). We address this challenge by proposing a structured semiparametric approach that helps modelers identify and calibrate an approximate set of low-dimensional deterministic calibration functions.
In this paper, we make a distinction between model variables (x) that can be controlled by the modeler (e.g., system geometry, and material choice) and the model parameters that are unique to the physical system being studied (e.g., Young's modulus, yield stress, and Poisson's ratio) [10]. In many cases, the true values of the model parameters are unknown as experiments are too expensive and/or the parameters cannot be measured directly. Calibration provides modelers with an approach to infer the values of these model parameters from two sets of data: (i) low-fidelity samples from a simulation model, and (ii) high-fidelity samples from either physical experiments or high-fidelity simulations [3,11]. This dataset can be used in a frequentist-based method that involves finding the values for the model parameters by minimizing the discrepancy between the low-fidelity and the high-fidelity models [11]. However, these methods generalize poorly as they do not provide confidence bounds on the identified model parameters. Alternatively, Bayesian approaches can address this limitation by quantifying the uncertainty in the model parameters [12–14].
Parametric functional calibration can be achieved by assuming the class of calibration functions a priori (i.e., ) and then calibrating the unknown calibration parameters [15]. The experimental efficiency of functional calibration can be improved using emulators and the sequential acquisition of training samples [16]. In addition, the authors of Ref. [17] considered a variant of functional calibration where the low-fidelity and the high-fidelity models do not share the same set of model parameters. Specifically, they proposed a linear input mapping approach on the inputs of the low-fidelity model and used a regularization term to limit over-fitting. This allows for flexibility in the class of functions but cannot be used for functional calibration when the inputs to the high-fidelity model are also inputs to the low-fidelity model. If the functional forms of the model parameters are not known, then a prior distribution can be placed on the space of unknown calibration functions , instead of the calibration parameters , and updating them using the available training data [18,19]. This approach is typically referred to as nonparametric functional calibration. In addition, the calibration accuracy of nonparametric calibration can be improved by formulating and solving it as a “shortest path” problem [20,21]. Despite these contributions, the computational cost of functional calibration is often too exhaustive for problems with a large number of calibration functions, nor does it provide support for what class of functions should be used.
In this work, we introduce a semiparametric calibration method to select an appropriate class of calibration functions aimed at reducing the number of calibration parameters and limiting the compromise in prediction accuracy. We are specifically interested in the functional calibration of a low-fidelity model to an intrinsically noisy high-fidelity data source. Intrinsic noise is a type of aleatoric uncertainty that cannot be reduced with the acquisition of additional data points and is encountered in simulation models (e.g., agent-based models) as well as in physical experiments [22,23]. For example, such problems arise during the calibration of low-fidelity CG molecular dynamics simulations, to high-fidelity AA models. After finding the calibration functions, the modeler can use the CG model to explore physical phenomena at lengths and time scales that are intractable with the AA model. Note that changing the lengths and time scales is an invasive modeling procedure that precludes the modeler from using emulation techniques and emphasizes the importance of having an appropriate deterministic approximation of the calibration functions. An earlier work that addresses this issue involves the use of a global optimization method coupled with a neural network to minimize the discrepancy between the CG and AA model by tuning the force field calibration parameters [24]. However, a drawback of this approach is that the modeler needs to define the functional form of the force field parameters a priori, risking poor calibration accuracy and overfitting the calibration functions. The Bayesian-based method presented in this work circumvents this issue through the use of uncertainty quantification (UQ) to inform what calibration functions can be replaced with low-dimensional approximations.
The semiparametric calibration method introduced in this work provides support in selecting an appropriate class of functions for problems with a large number of calibration functions (i.e., where ). We consider this to be a semiparametric approach as a satisfying class of calibration functions is selected within the method as opposed to defining it a piori. This is similar to the decision-support framework presented in Ref. [25] to identify and remedy interpretability issues encountered when using Bayesian calibration. First, we assume a high-dimensional and flexible class of candidate functions for each model parameter. This introduces a large number of calibration parameters for which we quantify their uncertainty through a modified approximate Bayesian computation scheme. Specifically, we use the Hessian of a distance function at its global optimum to reduce the number of samples necessary to infer the first two statistical moments of the calibration parameters. The uncertainty sources considered in this method include the intrinsic uncertainty in the computer simulations as well as the interpolation uncertainty associated with the use of emulators. Through the posterior distribution of the unknown calibration functions, the modeler gains insight into the sensitivity of each calibration function, their signal-to-noise ratio, and their correlation. This insight is then used to identify what calibration functions can be replaced with low-dimensional approximations, and what class of functions should be used. Consequently, we present a systematic approach that facilitates the synthesis of multiple sources of domain knowledge to support modeling decisions (e.g., data from physics-based simulations, experimental data, and expert knowledge).
In the remaining sections of this paper, we first introduce available parametric and nonparametric functional calibration methods (Sec. 2). Next, we introduce the newly developed semiparametric functional calibration method and its constructs (Sec. 3). To verify the method's performance and validity, we continue our discussion by testing it on a test problem and an engineering example, with three and 14 unknown calibration functions, respectively (Sec. 4). Finally, we conclude this work by summarizing its findings and introducing opportunities for future work (Sec. 5).
2 Background
Model calibration enables synthesis of data from: (i) simulation models with controllable degrees of fidelity [26–28], (ii) physical experiments and a single simulation model [12], and (iii) multiple simulation models [29]. In this section, we provide an introduction to the general formulation of available methods for model calibration and its extension to functional calibration.
2.1 Model Calibration.
where n is the number of observed samples, is a set of unknown but assumed to be constant modeling parameters, is a bias function that captures the systematic difference between the low-fidelity and the high-fidelity data source, ρ is a scaling parameter, and ei is a noise term. The modeling parameters in this formulation are typically associated with some physical constants but cannot be observed experimentally. In addition, note that the high-fidelity function has been written as a function of the modeling variable x. Technically, it is also a function of the modeling constants , but it is convention to not include this in the formulation as it cannot be observed experimentally, let alone controlled.
where ρ is assumed to be equal to 1, is the space of admissible values of the modeling parameters, and is the L2 norm. However, in the formulation of Eq. (2), the calibration parameters are assumed to be constant over the design space whereas many calibration problems show a functional relation with respect to the modeling variables [6,8,9].
2.2 Functional Model Calibration.
where is a normal distribution with mean and standard deviation . In addition, is a set of observations associated with a training dataset of modeling variables (and calibration parameters for the low-fidelity data source). The GP is typically modeled through a squared exponential or a Matérn kernel function [34]. Consequently, GP emulators provide a practical approach that can be used to replace costly low-fidelity and high-fidelity simulations to expedite the calibration process.
3 A Semi-Parametric Method for Calibration
In this section, we introduce a method for calibration function formulation, UQ and dimension reduction.
3.1 Overview of the Semi-Parametric Calibration Method.
The flow of the semiparametric calibration method is shown in Fig. 2 and consists of three steps: (i) data acquisition and emulation, (ii) UQ of candidate functions, (iii) and reduced dimensional calibration. The calibration method starts with the creation of a design of experiments of uniformly distributed input conditions for the low-fidelity and the high-fidelity model (Step 1a). A design of experiments with nl samples for the low-fidelity model is given as where is a d-dimensional vector of model variables and is a p-dimensional vector of model parameters. Note that this means that the experimental design of the low-fidelity model takes samples from the joint space . A design of experiments with nh samples for the high-fidelity model is given as where the vector is a d-dimensional vector of model variables defined on the same space as the low-fidelity model. Subsequently, we proceed by acquiring the response data for the low-fidelity model and for the high-fidelity model by performing the necessary simulations. To predict the response for unobserved combinations of model variables , and modeling parameter we train two GP emulators as and for the low-fidelity model and high-fidelity model, respectively (Step 1b). Subsequently, we use these emulators in an optimization scheme (Step 1c) to find the deterministic optimal calibration parameters for a candidate set with high-dimensional calibration functions that we assumed to be radial basis functions (RBFs) due to their flexibility and ease of implementation [14,18,31].
In contrast to finding a set of calibration parameters that minimizes the difference between the predicted mean responses, we want to infer the first two statistical moments of the calibration parameters as using a distance metric . This lends itself to an approximate Bayesian computation inference scheme that involves sampling the posterior distribution of calibration parameters (Step 2a) without the need for evaluating the likelihood function. This is a significant advantage as evaluating the likelihood function requires the inversion of the covariance matrix that quickly becomes intractable for problems with a moderate to large number of calibration functions (). In addition, we introduce a more data-efficient sampling process that uses the Hessian of the distance function to retrieve the eigenvectors of the calibration parameter covariance matrix (Step 2a and Step 2b). Subsequently, on each of these vectors, we find the two samples on the approximate Bayesian computation rejection boundary [35]. This enables a modeler to find two samples along each eigenvector that have the same distance function value (Step 2c). Subsequently, this provides a small yet efficient set of samples from which we approximate the first two statistical moments of the candidate set of calibration parameters' posterior distribution as and , respectively (Step 2d). This assumes that the distribution of the calibration parameter is jointly Gaussian, and as such, we used a hat to indicate that this is an approximation.
The posterior distribution of the calibration parameters can be used to draw insight into the functional form of the calibration functions. We achieve this by visualizing the posterior distribution of the calibration functions and by evaluating their signal-to-noise ratio (Step 3a). Subsequently, the confidence bounds, as well as the covariance matrix, provide insight into the underlying class of functions (Step 3b). If the confidence bounds are narrow and the mean function shows a strong trend (i.e., the expected range of the function is large) then this informs the modeler that the calibration functional forms need to be selected carefully. Moreover, in this scenario, the signal-to-noise ratio will be relatively large and a more complex calibration functional form (i.e., more input dimensions) can be selected. Conversely, if the confidence bounds are wide, or the trend is weak, then the signal-to-noise ratio will be small and the functional form of this calibration function is less consequential to the calibration accuracy and can thus be simplified (e.g., constant, or linear) without significantly affecting the calibration accuracy. Finally, the simplified set of calibration functions can then be used to find the optimal set of reduced dimensional calibration parameters . The presented method helps a modeler reduce the complexity of the calibration functions, provides insight into the physics that underlie their models, provide strong calibration accuracy, and enables the generalizability of these calibration functions to low-fidelity simulations at larger length and time scales. In the remainder of this section, we will provide additional details on the method's individual constructs: (i) distance and calibration functions, (ii) modified approximate Bayesian calibration scheme, and (iii) how to select low-dimensional calibration functions.
3.2 Formulation of Distance and Calibration Functions.
where is the range of the response variable. Note that in this expression we calculate the area under the curve of the probability predicted from the product of the low-fidelity emulator with the high-fidelity emulator. This has two desirable properties: (i) minimizing the distance function will force the mean predictions of the emulators to be the same, and (ii) the distance function places a higher importance on improving the goodness of fit in areas of the response surface with lower posterior predictive uncertainty. Another interesting observation is that the range of Eq. (5) is , where would correspond to a such that the two GP emulators are equivalent (i.e., and ). However, in most cases the data observed from the high-fidelity model is corrupted with noise (i.e., ) and as such the lower bound will be greater than zero. This can be corrected by taking the root of the term inside the integration of Eq. (5) and would provide a probabilistic interpretation for the distance function. However, doing so does not alter the calibration performance and thus we leave it as an interesting observation.
If we normalize the range of the model variables (i.e., ), then the aforementioned probabilistic interpretation of the distance function is retained.
where the ith element of the -dimensional vector k is defined as , the i, jth element of the matrix K is defined as , and c is an vector of weights (i.e., a set of parameters that need to be inferred from the training data). In this expression, the terms provide a set of nc radial functions centered at spatial locations . Taking the inner product of these radial functions with the weights c allows us to model a wide range of continuous functions. The nc center locations are best designed through a space-filling set of centers, in this work, we will use a uniform grid (i.e., a grid with a total of centers, with nc centers for each of the model variables that are an input to the ith calibration function). This introduces a set of calibration parameters to parameterize a total of calibration functions. The form of the calibration functions is equivalent to those used in Ref. [38], except that we use a smaller number of center locations. The big advantage of using RBF functions is that we can readily control the complexity of the calibration functions that we can capture by increasing the number of centers. In addition, using the RBF functions we can easily avoid extrapolating our emulation model by setting the range for the centers equivalent to the range of the associated dimension in the training dataset. In contrast, setting the right range for the constants of a set of polynomials is less straightforward [39].
3.3 Uncertainty Quantification of the Calibration Functions.
Using the distance metric in Eq. (6), we now want to directly sample the posterior distribution of the calibration functions through an approximate Bayesian computation scheme. This can be achieved through available methods but is inefficient for high-dimensional problems [40,41]. Consequently, we introduce an alternative approach that requires fewer samples of the posterior distribution of the calibration parameters. In Eq. (6), the GP emulators for the low-fidelity model and the high-fidelity model are used to create a set of pseudo-observations for any value of the calibration parameters [35]. Two typical approaches to sampling the posterior distribution of the calibration functions is the approximate Bayesian computation rejection algorithm [40] and the approximate Bayesian computation Markov chain Monte Carlo algorithm [41]. However, these approaches are known to be inefficient and intractable for problems with a large number of parameters [35,42]. In this subsection, we introduce a modified version of the approximate Bayesian computation rejection algorithm that reduces the number of necessary samples to approximate the first two statistical moments of the calibration function's posterior distribution. Specifically, the algorithm starts by identifying the eigenvectors of the posterior distribution of the calibration parameters and then finds the two samples over each vector that are at the rejection boundary.
The optimization of this function can be achieved through any gradient-based method starting from multiple initial calibration parameters or an evolutionary-based algorithm [43].
Subsequently, the eigenvectors can be obtained through an eigenvalue decomposition as , where the columns of Q are the eigenvectors , and is a matrix of mostly zeros except for its diagonals that contain the eigenvalues. Using the Hessian to infer the eigenvector is similar to taking the quadratic term of the local Taylor's series expansion at [43]. Moreover, it should be noted that using a Newton-type gradient-based optimization scheme (e.g., sequential quadratic programming) to find the global minimum of the distance function has the added advantage of providing an approximation of the Hessian H.
From the set of eigenvectors, we now need to infer the posterior distribution of the calibration functions (Step 2c and 2d). This can be done by finding the two samples with the same distance function values over each eigenvector. This idea has been visualized in Fig. 3 where the red lines indicate the eigenvectors that have been identified in Steps 2a and 2b and the surface plot represents the distance function given in Eq. (6). The algorithm then searches for two points over each eigenvector that have a distance function value equal to . In other words, we have selected the rejection boundary of the value to be half of the distance function response at the global optimum . This fraction can be changed but will mainly influence the magnitude of the covariance matrix. This is because the equal distance function boundary approximates an ellipsoid (as can be seen from Fig. 3). Note that each element of the sample set can be found by solving a one-dimensional optimization problem. Subsequently, we can approximate the first two statistical moments of the posterior distribution of the calibration parameters as and .
Note that this approach is a more data-efficient alternative to conventional approximate Bayesian computational rejection algorithm that would require significantly more samples from the prior distribution to approximate the posterior distribution. More specifically, by drawing samples from a noninformative prior distribution the conventional rejection algorithm would find a set of samples uniformly covering the elliptically shaped surface of the distance function in Fig. 3. Consequently, from this observation, we can conclude that the introduced sampling scheme will find a covariance matrix that would be the same as the approximate Bayesian computation rejection scheme with a uniform/uninformative prior. This implies that the introduced scheme does not enable a modeler to use an informative prior distribution for the calibration parameter, but we believe that this type of information will only rarely be available.
An analogy to the procedures in Steps 2a to 2d is that we are trying to find the minimum distance function value for which we can fit an -dimensional hyper-ellipse with a constant surface area over the distance function surface . This is an optimization problem that is remarkably similar to robust design optimization where modelers are concerned with finding a design that is robust with respect to various sources of uncertainty [44,45]. However, the difference with our approach is that the surface area of this hyper-ellipsoid is implicitly selected through the rejection boundary a priori.
3.4 Approximating the Class of Calibration Functions.
In the final steps of the method (Step 3a to 3c), we use the posterior predictive distribution of the unknown calibration functions to inform their parametric form. This process starts with propagating the uncertainty in the calibration parameters to quantify the uncertainty in the calibration functions . The quantified uncertainty and the range of the function can then be used to inform their parametrization. For example, consider the two plots in the bottom left of Fig. 2, here we have two functions with roughly the same magnitude of confidence bounds, but one shows a strong linear upward trend whereas the other appears to be constant. Consequently, the modeler has reason to believe that these functions can be approximated to be linear and constant, respectively. The simplified set of calibration functions can be established by considering multiple classes of low-dimensional functions (e.g., constant, linear, quadratic, sigmoid, and RBF). By comparing the relative difference between the uncertainty bounds and the shape of the candidate set of calibration functions the modeler can observe the most likely shape of each function and their relative importance to the calibration accuracy. However, visual inspection is not systematic and becomes complicated for calibration functions with more than two inputs (i.e., ).
where the numerator is the expected range of the ith unknown calibration function and the denominator is the maximum posterior standard deviation. A low signal-to-noise ratio (e.g., ) indicates poor identifiability as the range of the function is small compared to the range of the confidence bounds. Consequently, these functions can be simplified significantly as their form is of marginal consequence to the model's prediction accuracy. In contrast, if the signal-to-noise ratio is large (e.g., ), then this implies that the training data provides significant insight into the unknown calibration functions. Consequently, it is reasonable to select a more complex functional form that provides a good approximation of the unknown function. Finally, for calibration functions with a moderate signal to noise ratio (e.g., ) it is possible to glean some insight into their functional form. However, these functions are best described through low-dimensional approximations (e.g., linear or quadratic) as the “true” functional form is still likely obscured by the nontrivial confidence bounds. The distinction of the small, medium, and large signal-to-noise ratios has been established based on empirical results as presented in Sec. 4. A summary of recommendations to guide the selection of the appropriate functional forms is provided in Table 1. Alternatively, variance-based sensitivity analysis can be used to arrive at a metric to inform the class of calibration functions [46]. The advantage of variance-based sensitivity analysis is that it provides an approximation of the higher-order interaction effects. However, in this work we recommend the use of the signal-to-noise ratio as variance-based sensitivity analysis can be time intensive for high-dimensional problems and the interaction affects can also be gleaned from the off-diagonal element of the calibration parameter covariance matrix .
Signal to noise ratio | Small signal to noise | Medium signal to noise | Large signal to noise |
---|---|---|---|
Calibration functional form recommendation | Great potential for simplification, a constant function will be appropriate in most cases. | Some potential for simplification, a linear or quadratic function will be appropriate in most cases. | The calibration function will have to be selected carefully to match the semiparametric approximation. |
Additional considerations | This scenario indicates that the calibration functions are unidentifiable if the noise is relatively large. | These functions can be simplified with low-dimensional approximations and provide some insight into their “true” form. | The general trend of the approximate calibration function is likely consistent with the “true” form. |
Signal to noise ratio | Small signal to noise | Medium signal to noise | Large signal to noise |
---|---|---|---|
Calibration functional form recommendation | Great potential for simplification, a constant function will be appropriate in most cases. | Some potential for simplification, a linear or quadratic function will be appropriate in most cases. | The calibration function will have to be selected carefully to match the semiparametric approximation. |
Additional considerations | This scenario indicates that the calibration functions are unidentifiable if the noise is relatively large. | These functions can be simplified with low-dimensional approximations and provide some insight into their “true” form. | The general trend of the approximate calibration function is likely consistent with the “true” form. |
4 Results and Observations
The use of the introduced calibration method will be highlighted in this section by applying it to a test problem with three calibration functions, and the calibration of a low-fidelity CG simulation to a high-fidelity AA simulation with 14 unknown calibration functions.
4.1 Four-Dimensional Test Problem.
where . In order to make inferences into the form of the assumed to be unknown calibration functions , we start by creating a training dataset (Step 1a) including 100 samples of the high-fidelity model and 400 samples of the low-fidelity model. The samples of the low-fidelity model span a four-dimensional space that consists of the model variable x as well as the calibration parameters t. Subsequently, we proceed by training our emulation models and (Step 1b) and then use these to find the deterministic optimal set of parameters of the candidate set of calibration functions (i.e., using Eq. (6) as part of Step 1c). Moreover, we have set the number of RBF centers equal to five (i.e., nc = 5).
Using the approximate Bayesian computation scheme (Step 2) we acquire the posterior distribution of the calibration parameters and use it to get the posterior distribution of the calibration functions. Specifically, by drawing samples, we retrieve the calibration functions as presented by the purple shaded regions and the dashed-dotted lines in Figs. 4(a)–4(c). From these figures it can be observed that and have relatively small confidence bounds compared to . Moreover, this gives us the signal-to-noise ratios as , and . This informs us that selecting an appropriate functional form for and is more important than for . This observation is consistent with our expectation, as when we look at the high-fidelity function (Eq. (12)) we find that the third calibration function influences the frequency of a sinusoidal term and, as a result, variations in will decrease the value of Eq. (6) only marginally. More specifically, by observing the predicted response of the high-fidelity model and the low-fidelity model in Figs. 4(e) and 4(d), we find that the influence of the sinusoidal trend associated with the third calibration function is marginal. Moreover, if the frequency of this sinusoidal trend were to change as a result of , then, the confidence bounds will still overlap. However, changes in or are proportional to the deviation between the low-fidelity model and the high-fidelity model.
Comparing Fig. 4(a) with Fig. 4(b) we observe that the mean trend of the first calibration function (i.e., ) changes more significantly than the mean trend of the second calibration function (i.e., ). Subsequently, we have evidence that it is more reasonable to assume a more high-dimensional and complex function for the first calibration function (i.e., ). This is an encouraging observation as this is consistent with the relatively high signal-to-noise ratio and insights that we can draw from the problem formulation as stated in Eq. (12) (i.e., we know that the “true” functional form is logistic). Modeling such a function rigorously would require at least four calibration parameters. In contrast, shows a much simpler linear trend. Nevertheless, the confidence bounds can hide more complex functions, the only conclusion that we can make from this observation is that we can model the second calibration function with a linear class of functions and still expect reasonable calibration accuracy. This observation is even more relevant for the third calibration parameter where the confidence interval could hide a simple constant function or a wildly fluctuating multimodal function that would require many variables to parametrize. Nevertheless, these figures inform us that we can expect to have a reasonable calibration accuracy if we simplify the functional class of , and by an RBF with five centers a linear function and constant function, respectively. Moreover, note that these conclusions are consistent with the recommendations presented in Table 1.
With the simplified functional form of the three calibration functions , we can once again minimize Eq. (6) to get the reduced-dimensional set of optimal calibration parameters . This is a significant simplification in the number of parameters, as the RBF calibration with five centers for three calibration functions requires 18 parameters whereas the simplified set of functions only requires only nine parameters. In addition, the method allows a modeler to make an informed decision in terms of the class of functions to assume for each calibration function. If a modeler was to consider using a set of five functional classes (e.g., constant, linear, quadratic, cubic, RBF), then a problem with three calibration functions would have a total of alternative formulations. Exploring all these alternatives manually would be a tedious and time-consuming effort.
A concern that arises in the three-dimensional test problem is that the simplification will have a negative influence on calibration accuracy. When plotting the predicted response in Fig. 4(d) of the low-fidelity model and the high-fidelity model, we observe that it is difficult to separate the two from one another when using the high-dimensional functional parameterization (i.e., the dashed-dotted purple lines in Figs. 4(a)–4(c)). Contrasting this with that of the simplified functionalization (i.e., the solid black lines in Figs. 4(a)–4(c)) in Fig. 4(e), we find that the visual difference is hard to notice. In fact, the root mean squared percentage error (RMSPE) [47,48] for the original functionalization and the simplified functionalization are 0.36% and 0.38%, respectively. As might be expected, the calibration accuracy of the simplified functionalization is slightly worse, but still very good. In fact, an RMSPE less than 10% is typically regarded as excellent [49]. This is further corroborated by the observation that it is difficult to differentiate between Figs. 4(d) and 4(e).
In order to compare the tradeoff between the reduction in model complexity and calibration accuracy we can compare the results using the Bayesian information criterion that is given as , where is the likelihood and is the number of calibration parameters [50]. By calculating the Bayesian information criterion for observing the high-fidelity training dataset by sampling from the low-fidelity emulation model, we get BIC = 673.05 and for the mean approximation of the candidate set of functions and the simplified functions, respectively. The significance of this improvement should not be underestimated as it is on a logarithmic scale. Even though the Bayesian information criteria provide a statistically rigorous tradeoff between calibration accuracy and functional complexity, using it directly to optimize the objective would be complicated. This is because we would have to search over a design space with categorical alternatives, where nf is the number of considered functional classes.
4.2 Calibration of a Coarse-Grained Epoxy Model.
Modelers that are interested in exploring the thermomechanical properties of new epoxy-based composite materials frequently rely on simulation models to explore the high-dimensional design space of atomic structures and chemistries. The simulation models of choice are typically molecular dynamics; however, because AA simulations are intractable at experimentally verifiable length scales, modelers rely on CG simulations to predict thermomechanical properties. The force field of AA simulations is often well-defined, but CG models use empirical force fields with interaction parameters that must first be calibrated, for example, by matching the predictions of the higher-fidelity AA simulations. Moreover, previous studies have shown that the degree of crosslinking (DC) has a consequential effect on the properties of an epoxy system [51–56]. Despite this knowledge, previous efforts have been unable to discover the functional dependence of the model parameters on the DC. In this study, we will use the introduced calibration method to elucidate this functional dependence. The material science related findings of this study have already been published in a previous publication [57]. In contrast, in this paper, we demonstrate the design related findings.
We use an epoxy compound with Bisphenol A diglycidyl ether (DGEBA) as the resin and either 4,4-Diaminodicyclohexylmethane (PACM) or polyoxypropylene diamines (Jeffamine D400) as the curing agent. AA simulations are used as our high-fidelity model and CG simulations are used as our low-fidelity model. Moreover, we consider DC as the modeling variable, and each unique coarse-grained bead has two calibration functions, cohesive energy, and bead size. In the coarse grain procedure, we identified seven unique beads for both the PACM and D400 systems, and thus, we have a total of 14 calibration functions. We consider these calibration functions to be the inputs to the low-fidelity CG simulation that enables us to predict the modulus, yield strength, Debye-Waller factor , and the density of the two material systems (DGEBA + PACM and the DGEBA + D400). These properties have been selected as they are important markers of the overall thermomechanical properties. Using the introduced calibration method, we will simultaneously match all eight low-fidelity simulation model responses to the AA model responses. To calibrate the unknown functions with respect to multiple responses simultaneously we take their product at the inner integral of Eq. (6).
We started the calibration procedure by running 700 and 500 low-fidelity simulations for the DGEBA + PACM and DGEBA + D400 system, respectively. These simulations were performed using a fully sequential space-filling design obtain from a Sobol sequence [58]. The advantage of using a fully sequential space-filling design over a one-shot space-filling design is that we can sequentially add more simulation experiments until an emulation model with sufficient fidelity has been obtained [59,60]. In fact, this is the reason why 700 simulations where used to approximate the DGEBA + PACM CG model as opposed to the 500 used for the DGEBA + D400 system. In addition, we used a similar procedure for the high-fidelity model that included 19 simulations and 20 simulations for the DGEBA + PACM and DGEBA + D400 system, respectively. To quantify the uncertainty in the calibration parameters (Steps 2a to 2d) we assumed that for each of the 14 RBF functions we used three center points and one shared length-scale parameter. Subsequently, we have a total of calibration parameters for which we approximated their first two statistical moments. Using the approximate distribution of the calibration parameters, we quantified the uncertainty in the calibration functions as plotted in Fig. 5. In these plots, the purple shaded regions are the 95% confidence intervals and the purple dashed-dotted lines are the mean calibration functions. Moreover, the interaction energies are associated with the first seven calibration functions (i.e., ) and the bead sizes are associated with the last seven calibration functions (i.e., ). From these figures, we can infer that many of the calibration functions can be replaced with either a constant or linear functions. Specifically, using the signal-to-noise ration for each calibration function () we identified selected the functional form according to the followings conditions, if we used a constant function, for we have used a linear approximation, and for we have used an RBF with three centers. The selected and calibrated functions, as inferred from the training data, are shown by the black solid lines in Fig. 5. Finally, some simplified calibration functions display a nontrivial difference between the quantified confidence bounds. While this observation highlights the high-dimensional/complex interaction between the calibration functions, it also confirms that the simplified functions have a form that is consistent with the form of the all-RBF approximations.
The mean and the confidence interval for the eight CG model responses as predicted from the GP emulators for the simplified set of calibration parameters are plotted by the blue shaded region and lines in Fig. 6, respectively. The orange plots represent the high-fidelity simulation response. From these figures, we can observe that the response of the low-fidelity simulation model matches well with that of the high-fidelity response. In particular, we observe that the mean trend of the density appears to provide the best match to the AA model. The reason for this is that in the distance function in Eq. (6), we consider the predictive uncertainty of the GP emulators and as such the algorithm places more emphasis on those responses with smaller intrinsic modeling uncertainties.
In terms of fidelity, it would be of interest to compare the calibration accuracy of the simplified parametrization with the mean approximation of the all-RBF parametrization , and thus we have plotted the responses in Fig. 7. From these plots, we observe that the candidate set of all-RBF parametrization performs significantly better in terms of Young's modulus and the Yield stress. However, this comes at the expense of accuracy in the predicted Debye-Waller factor, . When comparing the accuracy of the mean prediction we find that the simplified parametrization and the candidate set of calibration functions have an RMSPE of 10.0% and 9.8%, respectively. As expected, the simplified parameterization performs slightly worse, but this change is only marginal compared to the decrease in the number of parameters necessary to parametrize the calibration functions.
The simplified set of calibration functions contains 1 RBF, 4 linear, and 9 constant functions for a total of 21 calibration parameters. This is a significant reduction compared to the 43 parameters from which we started. In addition, if we were to consider that each calibration function could be modeled through one of five functions then the modeler would have to choose one class of functions from a set of 514 different calibration function combinations. Exploring more than 6 billion alternative sets of calibration functions is too time-consuming and as such the introduced calibration method provides a modeler with valuable support. The simplified set of calibration functions have been selected to reduce the number of calibration parameters and avoid over/under-fitting so that they can be used to reliable explore the properties of larger material systems using the low-fidelity CG model.
A final insight that can be obtained from this analysis is that the average absolute correlation between the calibration parameters of different calibration functions as plotted in Fig. 8. The conclusion that can be drawn from this figure is that and have a high correlation with many other calibration functions. This observation is in line with what is expected as these calibration functions are associated with the cohesive energy of the coarse-grained beads that are found in both the DGEBA + PACM and the DGEBA + D400 system. In addition, we can observe that the correlation between the calibration functions associated with the size of the beads (i.e., ) are correlated among themselves but have little correlation with the interaction energies. This should come as little surprise as the density mainly depends on the bead size and the other mechanical properties depend on the cohesive energy. We make this observation here as it validates the potential of the introduced calibration method, for a more detailed discussion on the modeling aspects of this study we refer the reader to a paired study [57].
4.3 Validation of Selected Calibration Functional Form.
A limitation of the introduced semiparametric calibration scheme is that it requires the modeler to select the appropriate calibration functional form based on the signal-to-noise ratios. As a consequence, given the same training dataset, different modelers can arrive at different sets of calibration functions . We argue that this is an advantage of the proposed method, as it enables the modeler to leverage their implicit knowledge when selecting the calibration functional form. However, a counterargument would be that in the ideal scenario, this knowledge is explicitly accounted for by the physics of the simulation models and that the calibration functions are selected autonomously. Consequently, to provide additional corroboration of the introduced semiparametric calibration method, it would be beneficial to compare it to an approach where the selection of the calibration functional form is automated.
where the calibration functions are given as and is the ith response of the jth high-fidelity observations. Note that in Eq. (13) the regularization is only with respect to the parameters associated with the linear terms of the polynomials (i.e., ). Consequently, by minimizing Eq. (13) we only penalize functions that deviate from being constant. The parameters are regularized with respect to the L2 norm, the advantage of this is that we have a continuous objective function. However, the disadvantage is that in the optimization scheme parameters that are close to zero have a small gradient, and thus might not converge to zero as they would using the discontinuous L1 norm. Finally, we believe that this is an apt approach for validation as minimizing Eq. (13) can be expected to be on par with the introduced semiparametric calibration approach. The reason for this is that in the final set of calibration functions only one was not linear or constant. While the degree of model simplification is problem dependent, the primary advantage of the proposed method is that it informs when simplification of the calibration functions is realistic (i.e., when the signal-to-noise ratio is small).
A feature of Eq. (13) is that we can control the tradeoff between calibration accuracy and model sparsity through the regularization term γ. Setting a small value for γ will minimize the penalty on nonconstant calibration functions and thus provide more flexibility. In contrast, selecting a large value for γ will emphasize model sparsity by penalizing nonconstant calibration functions. To investigate the influence of the regularization term, we have calibrated the CG model using Eq. (13) for uniformly spaced values of . The resulting calibration accuracy has been plotted in Fig. 9(a) and the model complexity has been given in Fig. 9(b). The model complexity is the sum of the number of constant terms and the number of linear terms with a value greater than 1% of the admissible range of the calibration functions. Specifically, is an identity function where is the ith normalized nonlinear parameter. The discrete calibration results for 100 values of γ are presented by the blue circles and the blue lines present a moving average with a moving window of five. For comparison we also included the results obtained from the introduced semiparametric calibration approach, see the red dashed line. In addition, the black dotted line represents the first set of calibration functions obtained from Eq. (13) that was equal in model complexity to the results obtained from the introduced semiparametric calibration approach.
It is encouraging to observe that only considering constant and nonlinear calibration functions nearly doubles the RMSPE compared to the introduced semiparametric calibration method. Moreover, it shows that considering more complex functions is a valuable exercise as it can greatly improve calibration accuracy without compromising model complexity. Specifically, from these observations we can make two conclusions:
Decreasing γ increases model complexity and results in improved calibration accuracy. However, the calibration accuracy is not reduced to zero and as such considering only constant and linear functions does not provide sufficient flexibility.
Increasing γ decreases model complexity and calibration accuracy. However, the model complexity does not reduce to zero, and as such it would be valuable to use an L1 norm in Eq. (13) to truly enable the selection of the functional form.
Nevertheless, the introduced semiparametric calibration method provides a valuable method to balance model complexity with calibration accuracy albeit dependent on the modeler's intuition or expert knowledge.
From the previous study, we can also learn the importance of the various calibration functions and their form. Specifically, we plotted the normalized linear terms (i.e., ) of the Interaction energies and the bead size with respect to different values of the regularization term γ in Figs. 10(a) and 10(b), respectively. What we observe from these figures is that the parameters with larger values over the entire range of γ are the more important parameters, and functions, to achieve high calibration accuracy. Specifically, we find that ϵ2, ϵ3, and ϵ4 are highly influential to calibration accuracy. This is consistent with our observations in Fig. 5 where ϵ2, ϵ3 are approximated through a linear and RBF function, respectively. While ϵ4 has been approximated through a constant function, something that can be explained by noting that ϵ4 is highly correlated with ϵ3, as observed in Fig. 8, and its relatively small signal-to-noise ratio. Finally, we find that it is more important for the interaction energies to have more flexibility compared to the bead size. In summary, we find that we better calibration results can be achieved with the introduced semiparametric calibration method and also provides additional insight into the unknown calibration functions.
While using Eq. (13) provides validation for the introduced semiparametric calibration approach, it also stipulates a path forward. Calibration results can be improved by providing a flexible class of calibration functions in Eq. (13) and by enabling autonomous selection. However, doing so is a nontrivial challenge, and as such, this will be of interest if future work.
5 Conclusions
In this paper, we presented a semiparametric functional calibration method for models with input-dependent parameters. Specifically, we introduced a modified Bayesian approximation scheme that uses the Hessian of a distance function to quantify the uncertainty in a high-dimensional candidate set of calibration functions. Subsequently, the quantified uncertainty was then used to define a more appropriate class of calibration functions with fewer parameters. The novelty of the presented method is the joint consideration of the following three modeling characteristics: (i) it allows the modeler to benefit from various sources of domain knowledge (e.g., the physics embedded in their simulations models, data from physical experiments, and expert knowledge), (ii) the modeler can greatly reduce the complexity of the calibration functions while maintaining strong predictive accuracy for high-dimensional problems (i.e., more than 10 calibration functions), and (iii) it can be considered as a semiparametric approach as it does not require the modeler to define the form of the calibration functions a priori. These characteristics make the introduced calibration method a powerful tool for modelers to explore physical phenomena at lengths and time scales unattainable with their high-fidelity data sources. This could have broad reaching implications on materials design, environmental modeling, and operations engineering. Finally, we have demonstrated the validity of the proposed method through an engineering problem where the number of calibration parameters was reduced by over 50% while only increasing the predictive error by 0.2%.
Future research in nonparametric functional calibration can be directed toward the consideration of a systematic bias between the low-fidelity and the high-fidelity data sources. Simulations often do not account for all the physics associated with the physical processes that they aim to model. Consequently, a systematic bias exists between the low-fidelity simulations and the high-fidelity physical experiments, necessitating the consideration of a bias function. However, the bias function introduces significant modeling flexibility that is likely to lead to identifiability issues. Consequently, a second direction for future research includes investigation into the calibration function identifiability. Available methods that address the identifiability issue involve stringent assumptions on the form of the bias function, necessitate an informative prior distribution, or require additional response variables that might not always be available. Third, the generalizability of the developed method could be tested by addressing challenges in other disciplines. For example, the modified approximate Bayesian computation scheme could be used to quantify the prediction uncertainty in the hyperparameters of deep learning models. Further research in Bayesian-based functional calibration will enable modelers to get a de-eper understanding of the systems that they are studying, and thus facilitate the design of products with improved utility.
Acknowledgment
Support from AFOSRFA9550-18-1-0381, U.S. Department of Commerce under Award No. 70 NANB19H005 and National Institute of Standards and Technology as part of the Center for Hierarchical Materials Design (CHiMaD). In addition, versions of Figs. 5–8 have originally been published in NPJ Computational Materials [57] under the creative commons attribution-noncommercial-sharealike license.
Funding Data
Air Force Office of Scientific Research (Award No. 70 NANB19H005; Funder ID: 10.13039/100000181).
National Institute of Standards and Technology (Award No. Chimad; Funder ID: 10.13039/100000161).
Data Availability Statement
The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.
Nomenclature
- =
distance metric function
- d, p =
number of modeling variables and number of modeling parameters
- =
expectation and covariance function
- f =
function and function value
- =
Vector of kernel function values, matrix of kernel function values, length scale parameter, and center values of RBF functions
- =
Eigenvector and column matrix of eigenvectors
- n =
number of samples
- =
normal distribution
- P =
probability
- =
a matrix of left singular vectors, a matrix of right singular vectors, and the matrix of singular values
- X =
set of samples at the approximate Bayesian computation rejection boundary
- γ =
smoothing parameter for regularization
- ρ =
scaling parameter in the general calibration formulation
- δ =
bias function
- =
covariance matrix
- =
space of admissible calibration functions
- =
space of possible function responses
- =
set of calibration functions
- =
set of calibration parameters
- =
reduced set of calibration parameters
- ξ =
approximate Bayesian computation rejection value
- =
mean and standard deviation
- =
space of admissible calibration parameters and reduced dimensional calibration parameters
- =
space of admissible modeling variables and modeling parameters
- =
noise terms
- (k) =
index for the kth eigenvector
- =
number of updates and total number of update of the sequential optimization scheme
- =
simplified approximation
- –1 =
inverse of a square matrix
- =
an optimized parameter or function
- T =
vector transpose
- =
approximated parameters or function
- =
a variable arbitrary fixed at a specific value
Appendix: Discussion on the Introduced Method and Its Implementation
Concerning the proposed method's theoretical validity, as defined in Ref. [61], we continue our discussion by exploring some of its characteristics, limitations, and user scenarios in this section.
In the introduced calibration method we approximate the first two statistical moments of the calibration functions. The uncertainty in these functions comes from the intrinsic data uncertainties as well as the interpolation uncertainty. By using the Hessian of the distance function at the global optimum we infer the covariance matrix of the calibration function parameters. However, in this process, we search for a rejection boundary that has a distance function value (Step 2c). This is a challenge also observed in approximate Bayesian computation [35] and only allows us to find a scaled version of the covariance matrix and as such also scales the range of the confidence bound of the approximated calibration functions. Consequently, the calibration method will only provide insight into the relative importance of calibration functions as opposed to absolute importance.
An additional consideration is the potential multi-modality of the distance function. The introduced uncertainty quantification scheme infers the rejection boundary of the response around the global optimum of the calibration parameters and as such is unable to consider the local optimums isolated from the global optimum. Note that we used the term isolated, as the local optimum that is in the proximity of the global optimum will be explored by the calibration method. However, we believe that the influence of this limitation is marginal as the local minimum observed in our test problems has relatively small distance function values compared to the global optimum and as such would have a very small probability in a full Bayesian scheme (i.e., non approximate).
The introduced method is most suitable for the scenario where we have a low-fidelity model and a high-fidelity model that given the true set of calibration functions has no systematic bias. When used for two simulation models this can be a reasonable assumption as they can be constructed to include the same set of physics. However, when we extend this to the calibration of a low-fidelity simulation model to a set of high-fidelity physical experiments, then the existence of a bias function needs to be considered. The reason for this is that those simulation models typically include only the most prevalent physics and as such are expected to have a systematic bias with physical experiments [12]. Modelers need to take care when using the introduced calibration method in such scenarios as the analysis will try to compensate for the missing physics by using the additional freedom provided by the candidate set of calibration functions.
Another important consideration is that the calibration method requires the modeler to define a range for the calibration functions a priori. If this range is chosen to be too small then the true calibration function might extrapolate from this range and as such the identified calibration will provide only a little, or even misleading, information about the “true” form of the calibration functions. However, in most cases, this scenario is identified by finding a set of calibration functions that are at the calibration function bounds. In such a scenario, the modeler is advised to run additional low-fidelity simulations to increase the admissible range of the calibration functions.
A final issue to discuss is identifiability [62,63], a phenomenon that occurs when there is either a non-unique set of optimal calibration functions or when there is a large range for the calibration function parameters that have similar calibration accuracy. The issue of identifiability can be observed in the introduced method when the variance (i.e., the diagonal elements of the calibration parameter matrix) is significantly larger than that of the other parameters (e.g., compared to η13 in the coarse-grained epoxy model as studied in Sec. 4.2). In such an event, the associate calibration function can in most cases be replaced by a constant function or linear function (e.g., and in Fig. 5). However, this does not imply that the “true” calibration function has a low-dimensional parametric form, we can only conclude that its form has little influence on the overall calibration accuracy compared to other calibration functions. We have found in the coarse-grained epoxy example that the identifiability of the model improved when calibrating all eight responses simultaneously as opposed to individually. This is an observation that makes intuitive sense as including multiple responses reduces the likelihood that there exists a non-unique solution [64].
The purpose of discussing the above issues is to bring attention to the various assumptions and considerations that underpin the presented calibration method. Moreover, knowledge of these issues informs how results should be interpreted and guides further scientific effort.