Information fusion and machine learning for sensitivity analysis using physics knowledge and experimental data Thanks: Citation: Kapusuzoglu, B., & Mahadevan, S. Information fusion and machine learning for sensitivity analysis using physics knowledge and experimental data. Reliability Engineering & System Safety 214, 107712 (2021). DOI: 10.1016/j.ress.2021.107712

A direct strategy to enforce physics constraints in ML model predictions is by including the constraints within the loss function used in training the ML model [26]. Thus, while training a PIML model with inputs $\mathbf{X}_{obs}$ and outputs $\mathbf{Y}_{obs}$, the physical constraints can be incorporated as additional penalty terms in the loss function:

where $\mathcal{L}_{\rm ML}$ is the log marginal likelihood of the data for a GP model:

with experimental observations {$\mathbf{X}_{obs},\mathbf{Y}_{obs}$} being the training data for the GP model. For a DNN $\mathcal{L}_{\rm ML}$ represents training a loss function that evaluates a supervised error, e.g., root mean squared error (RMSE):

which measures the accuracy of predictions $\mathbf{\hat{Y}}$ for $n$ training samples. (Note that for the GP model, the likelihood is maximized, whereas for the DNN model, the RMSE is minimized). The additional physics constraint loss function $\mathcal{L}_{\rm phy}$ in the second term of Eq. (9) is weighted by a hyperparameter $\lambda_{\rm phy}$; the value of $\lambda_{\rm phy}$ controls the strength of the physics constraint enforcement. The inclusion of the second term ensures physically consistent model predictions and helps to reduce the generalization error, which is a measure of how accurately a model is able to predict the output QoI for previously unseen data [26].

The physical inconsistencies in the model predictions are evaluated using the physics constraint loss term. The generic forms of these physical relationships can be expressed using the following constraints:

where $\mathbf{\Gamma}$ denotes other variables or thresholds that define the physics constraint regarding the model output $\mathbf{\hat{Y}}$. Equation. 12 indicates that the constraints may take the form of equalities or inequalities. These equations can involve algebraic relationships or partial differentials of $\mathbf{\hat{Y}}$ and/or $\mathbf{\Gamma}$. The physics-based loss functions for these equations can be defined as:

where $\rm{ReLU}(x)=max(0,x)$ represents the rectified linear unit function and it acquires value when a threshold is violated in the inequality constraint, i.e., it penalizes the optimization when $\mathcal{F}_{2}>0$. It can also be used to penalize deviations from a desired physically consistent relationship among multiple outputs $\mathbf{\hat{Y}}$ [38].

The ML model output accuracy and uncertainty are dependent on the quality and quantity of the available training data. In some systems, the high cost associated with conducting experiments makes it infeasible to have adequate amount of training data to build purely data-driven models. Thus, it may be desirable to combine the physics-based model and available experimental data in seeking to maximize the accuracy of the sensitivity estimates. When the experiments are expensive, they can only be conducted for a few values of the inputs, whereas it might be possible to run the physics-based model for a larger set of input values. In that case, the simulation data can be used to pre-train an ML model, which is used as the initial model to be updated with experimental observations. Further, training of ML models requires the choice of initial values of the model parameters. The transfer of physical knowledge using a pre-trained ML model can prevent poor initialization due to lack of knowledge regarding the initial choice of ML model parameters prior to training.

Since the pre-training based on the physics-based model can use a large amount of training data (with multiple input parameter combinations) over a wide range of values, the pre-training may also help the eventual ML model to have wider generalization beyond experimental data. In the numerical example in Section 4, the pre-training strategy exercises the physics model over 1310 input combinations, whereas only 39 experiments are available. However, if the physics model is computationally expensive, then the advantage of the pre-training strategy in using a larger input data set (for physics model runs) compared to the experiments becomes limited.

The two proposed strategies to predict the QoI are shown in Fig. 2. Figure 2(a) shows the first method, where the physical knowledge is included through constraints within the loss function of an ML trained with only experimental data. Figure 2(b) shows the second method, where an ML model is first trained with data generated using the physics-based model and then updated using experimental data. Figures 2(c) and 2(d) show the trained ML model predictions ($\boldsymbol{\hat{Y}}$) for the two proposed strategies, respectively. The proposed PIML strategies can be applied to any physical system by leveraging the related physical constraints or physics-based models.

In Model 1, denoted as $\mathbf{\rm GP}$, only experimental observations are used for training. The hyperparameters of the GP model (process variance, correlation length scale along each input dimension, and trend function coefficients, and also measurement error variance if unknown) are optimized during training by maximizing the log marginal likelihood function shown in Eq. (8). In calculating the likelihood, the difference between the true response of the system $\mathbf{Y}_{\rm true}$ and the observed response $\mathbf{Y_{\mathrm{obs}}}$ is attributed to the observation error $\epsilon_{\rm obs}$, which is often treated as a zero-mean Gaussian random variable with variance $\sigma_{\rm obs}^{2}$.

Model 2, denoted as $\mathbf{\rm GP^{\mathcal{L}_{\rm phy}}}$, incorporates the first PIML strategy by enforcing physics constraints during the optimization of the GP model hyperparameters. More specifically, the physics constraints are included during the maximization of the log marginal likelihood function (Eq. (8)) while inferring the hyperparameters of the GP model. Thus, the training of Model 2 is achieved by maximizing the function in Eq. (14):

where $\mathbf{\hat{Y}}$ is the GP model prediction. Note that since $\mathcal{L}_{\rm GP}$ is to be maximized, the second term corresponding to the physics constraint has a negative sign. Gaussian process modeling under constraints has been studied in the literature [8, 9, 15, 32, 44]. Veiga et al.  [9] developed a framework that incorporates bound, monotonicity and convexity constraints in GP modeling. Golchi et al. [15] developed a Bayesian approach to GP modeling that incorporates the monotonicity constraint. The need to obtain the monotonicity information at each of the points in the derivative input set can slow down the computation as the input dimension increases since the size of the covariance matrix depends on the input dimension. This paper does not use the methods described above. Instead, the paper proposes a different method that penalizes violations of the physics constraints by introducing a regularization term in the likelihood function. To the best of our knowledge, this is the first study to apply the penalty approach to the likelihood function of the GP model. Further, the computational effort of the proposed method (i.e., the calculation of the regularization term) does not increase with the problem size since the regularization term does not need the inverse of the covariance matrix.

The current work considers two different approaches for the second PIML strategy. Both approaches use a pre-trained model, obtained using the physics model input-output data. In the first approach, the model parameters are updated using the experimental data; and in the second approach, a discrepancy correction term is added to the pre-trained model. The first step of pre-training using data generated by the physics model can be thought of as similar to a lower fidelity model, and the second step of improving the pre-trained model (either by parameter updating or by adding a discrepancy term) can be thought of as similar to incorporating higher fidelity data (experimental data, in this case) to improve the model. Various multi-fidelity modeling strategies with different combinations of the low- and high-fidelity models have been studied in the literature, such as filtering, fusion, and adaptation [42]. The discrepancy correction approach pursued here adopts the simplest strategy, namely, additive correction, where a model discrepancy term is added to the low-fidelity model [1]:

where $f_{HF}(\cdot)$ and $f_{LF}(\cdot)$ are the high and low fidelity models, respectively, $\mathbf{X}$ is the input to the model, and $\theta_{\delta}$ are the parameters of the discrepancy correction term. The correction term can be obtained using any suitable surrogate modeling technique.

Model 3, denoted as $\mathbf{\rm GP^{\rm upd}}$, pursues the first approach of the second PIML strategy, where a GP model is pre-trained using the coupled multi-physics model input-output and then updated with experimental data. Then, the model parameters of this pre-trained network are updated using the experimental data.

An alternative approach, denoted as $\mathbf{\rm GP^{\rm MF}}$, pursues the second approach of the second PIML strategy, i.e., it pre-trains a GP surrogate model with data generated from the physics model, then improves the surrogate using experimental data. Consider a physics model $G(\boldsymbol{\cdot})$ that maps input variables $\mathbf{X}$ and model parameters $\boldsymbol{\theta_{m}}$ to the numerical model output $\mathbf{Y}_{m}$:

Let $n_{D}$ be the number of collected observation data $\mathbf{Y_{\mathrm{obs}}}$ from experiments with input variable settings $\mathbf{x}^{(1)},...,\mathbf{x}^{(n_{D})}$, where $\mathbf{x}^{(i)}$ is the input variable setting for the $i$th experiment. The physics model prediction is inaccurate due to missing physics or due to other approximations. Thus, a model discrepancy term $\boldsymbol{\delta}(\mathbf{X})$ as a function of model inputs is introduced to capture the difference between $\mathbf{Y}_{\rm true}$ and $\mathbf{Y}_{m}$ [27]:

The experimental observations $\mathbf{Y}_{\rm obs}$ can be described in terms of the true system response $\mathbf{Y}_{\rm true}$ and the corresponding observation errors $\epsilon_{\rm obs}(\mathbf{X})$ as

Combining Eqs. (17) and (18),

When the physics model is computationally expensive, it is replaced by a cheaper surrogate model. In Model 3, a GP surrogate model is used to approximate the original physics model. The accuracy of the surrogate model prediction depends on the quality and quantity of the training data generated by the original physics model. The surrogate model error ($\epsilon_{\delta}(\mathbf{X})$) can be incorporated as follows:

where $\mathbf{\hat{Y}}_{m}$ is the surrogate model prediction.

A common approach to estimate the discrepancy term $\boldsymbol{\delta}(\mathbf{X})$ is the one formulated by Kennedy and O’Hagan [27], which is applicable in the context of Bayesian calibration. In that case, physics model parameters are sought to be calibrated, and a discrepancy term is added in the calibration equation. The discrepancy term can be expressed in multiple ways, such as constant, Gaussian random variable with unknown parameters (either input-dependent or not), or Gaussian process (either stationary or non-stationary) [31]. The hyperparameters of the discrepancy term are then estimated along with the physics model parameters using Bayesian calibration [25].

However, the situation considered here is much simpler. There is no calibration of the physics model parameters here; only the discrepancy term is needed. (In other words, the physics model parameters are already established). In that case, the model discrepancy can be evaluated for different input values of experimental tests and realizations of observation errors as follows:

Moving the surrogate model error $\epsilon_{\delta}(\mathbf{X})$ to the left-hand side, we can express the difference between the actual response and GP model prediction as

In this work, a second GP model is trained for $\boldsymbol{\hat{\delta}}(\mathbf{X})$ in terms of the inputs. Thus two GP models are trained in Model 3. The first GP model is constructed using the physics model input-output data, and predicts $\mathbf{\hat{Y}}_{m}$. The second GP model is constructed using the experimental data and the corresponding surrogate model predictions, and predicts $\boldsymbol{\hat{\delta}}$ (the difference between the surrogate model prediction and actual system response).

The GP model for the model discrepancy captures the combined contribution of measurement error, physics and surrogate model errors for a given prediction. Thus, the predictions of the first GP model (pre-trained) are corrected with the second GP model predictions ($\boldsymbol{\hat{\delta}}$) representing the model discrepancy term and can be written as

Model 4, denoted as $\mathbf{\rm GP^{\rm upd,\mathcal{L}_{\rm phy}}}$, combines both PIML strategies for GP, where the optimized model parameters of the pre-trained model based on the physics model input-output are used as the initial values. These model parameters are updated using the experimental data by minimizing the augmented loss function shown in Eq. (14) which consists of both the training loss function and the physics constraint loss terms.

In Model 5, denoted simply as $\mathbf{\rm DNN}$, a deep neural network is trained using only experimental data. In order to train the model, an optimization algorithm is used to find a set of model parameters (weights and biases) that best map inputs to outputs. The number of epochs, which is the number of complete passes through a batch of training dataset, layers and neurons needs to be optimized to reduce cost and improve model accuracy. The model accuracy is evaluated using a validation dataset not used for training; the number of epochs is gradually increased until the accuracy improvement is insignificant. The accuracy and generalization performance of the DNN are also affected by the dropout rate used in training the model; the dropout concept will be discussed in Section 3.3 below, and the selection of the dropout rate will be discussed in Section 4.

Model 6, denoted as $\mathbf{\rm DNN^{\mathcal{L}_{\rm phy}}}$, extends Model 5 by implementing the first PIML strategy, i.e., physical knowledge related to the physical process is enforced through constraints within the loss function of the DNN, as shown in Eq. (9). The physics-based loss function terms are evaluated for given experimental inputs at every optimization iteration during the training of $\mathbf{\rm DNN^{\mathcal{L}_{\rm phy}}}$; this makes the optimization process slower during the training. In particular, the multipliers in the physics constraint penalty terms ($\lambda^{\rm DNN}_{\rm phy}$) affect the training speed; a higher value of the multiplier makes the penalty stronger (i.e., the physics constraints more stringent), thus requiring an increased number of iterations to converge to the optimum solution.

Model 7, denoted as $\mathbf{\rm DNN^{\rm upd}}$, pursues the second PIML strategy, where a DNN model is pre-trained using the coupled multi-physics model input-output and then updated with experimental data. The pre-trained model is first trained with physics model input-output data consisting of input combinations over a range of values. Then, the weights and biases (model parameters) of this pre-trained network are updated using the experimental data.

Model 8, denoted as $\mathbf{\rm DNN^{\rm upd,\mathcal{L}_{\rm phy}}}$, combines both PIML strategies for DNN, where the optimized model parameters of the pre-trained model based on the physics model input-output are used as the initial values. These model parameters are updated using the experimental data by minimizing the augmented loss function shown in Eq. (11) which consists of both the training loss function and the physics constraint loss terms. Similar to Model 6, the inclusion of physics constraints makes it slower for the optimization to converge to optimal model parameter values.

An additive manufacturing application is used to illustrate the proposed PIML models for GSA and compare their performance. A fused filament fabrication (FFF) process is considered; commercial material Ultimaker Black Acrylonitrile butadiene styrene (ABS) is extruded from an Ultimaker 2 Extended+ 3D printer to manufacture parts with unidirectionally aligned filaments, and the porosity of the manufactured part is measured. FFF is a widely used additive manufacturing (AM) process due to its easy operation, low cost, and suitability for complex geometries. As the molten filament is deposited layer upon layer through a nozzle, it cools down, solidifies and bonds with the adjacent filaments. Rectangular specimens of length 35 mm, width 12 mm, and thickness 4.2 mm are manufactured for the ABS amorphous polymer, with constant filament height, width and length (0.7, 0.8, and 35 mm, respectively).

The output QoI is the porosity of the printed part, and the inputs are two process parameters, namely nozzle temperature and speed. The porosity of an FFF part is dependent on the temperature history at the interfaces between filaments. Thus, it is important to predict the temperature evolution of filaments for estimating the final mesostructure of the printed part. The analytical solution proposed by Costa et al. [7] for transient heat transfer during the printing process in FFF is used to predict the temperature evolution of filaments. A physics-based sintering model is developed, which considers realistic filament geometry, and allows the filament geometry to change during the printing process [24]. This model is used to predict the porosity of the FFF part using the temperature evolution of filaments, material properties, part geometry, and process parameters as inputs. Thus the mapping from input to output is a multi-physics model, i.e., models of two physical phenomena (heat transfer and sintering) are combined to predict the porosity given the values of two process parameters: nozzle temperature and nozzle speed.

The statistical properties of the QoI were observed to have negligible variability along the length of the specimens; therefore only the porosity measurements taken at the midpoint cross-section (see Fig. 3) are discussed here. These measurements were based on microscopy images processed through the ImageJ software [49]. Filaments were extruded through a nozzle with 0.8 mm diameter. The build plate temperature was constant and set to $110^{\circ}$C. Using Latin hypercube sampling, 39 sets of process parameters $\mathbf{X}$ were generated, and experiments were conducted at these 39 values. The ranges considered for the two process variables were: nozzle temperature $T$: (210${}^{\circ}C$ - $260^{\circ}C$), and nozzle speed $S$: (15 mm/s - 46 mm/s).

The basic ML models, namely Model 1 (for GP) and Model 5 (for DNN) are simply trained with the 39 sets of process inputs (temperature and speed) and output (porosity).

In the training of Model 2 ($\mathbf{\rm GP^{\mathcal{L}_{\rm phy}}}$) and Model 6 ($\mathbf{\rm DNN^{\mathcal{L}_{\rm phy}}}$), we impose two physics constraints (i.e., two separate loss function terms, $\mathcal{L}_{{\rm phy},k}(\mathbf{\hat{Y}})$, where $k=\{1,2\}$ and $\mathbf{\hat{Y}}$ is the porosity prediction). The corresponding loss function terms are defined as

considering physics violations related to the porosity in all the $N$ samples. In the first loss function, a negative value of porosity is treated as a physics violation. The second loss function penalizes the model when the predicted final porosity $\mathnormal{\hat{Y}}_{i}$ is greater than the initial porosity $\phi_{0,i}$ of the $i$th part. This is based on the physics knowledge that the total void area decreases as the bond formation takes place. Thus, the porosity predictions are ensured to be physically meaningful with the inclusion of these physics-based penalty terms.

The overall “loss” function of the GP model is

Note that the GP model parameters are obtained by maximizing the above function.

The overall loss function of the DNN model is

Note that the DNN model parameters are obtained by minimizing the above function.

In Model 3 ($\mathbf{\rm GP^{\rm upd}}$) and Model 7 ($\mathbf{\rm DNN^{\rm upd}}$), the ML models are pre-trained using the multi-physics model input-output. The pre-trained ML models are then updated using the experimental data. The training data for pre-training consists of 1310 input parameter combinations over a range of experimental values, i.e., (210${}^{\circ}C$ $\leq T\leq$ $260^{\circ}C$, 15 mm/s $\leq S\leq$ 46 mm/s ). Note that there are only 39 physical experiments available; this is one of the advantages of the pre-training/updating strategy, where the pre-training can be over a much larger set of input combinations, thus improving the generalization performance of the updated model. The input data are normalized prior to the training of the ML models (i.e., the output quantity porosity is dimensionless and between 0 and 1).

Model 4 ($\mathbf{\rm GP^{\rm upd,\mathcal{L}_{\rm phy}}}$) combines both PIML strategies for GP, and consists of two GP models: (i) the first GP model is trained using the physics model input-output samples consisting of 1310 input parameter combinations; and (ii) the second GP model is built for the discrepancy between the first GP model prediction and the actual system response using the experimental data, by maximizing the function shown in Eq. (35) to optimize the hyperparameters of the second GP model.

Model 8 $\mathbf{\rm DNN^{\rm upd,\mathcal{L}_{\rm phy}}}$, which is a combination of the two PIML strategies, uses the DNN model parameters trained using the physics model input-output as the initial values. Then, during the updating phase with the experimental data, these parameters are updated by minimizing the loss function shown in Eq. (36).

The four GP models (Models 1 to 4) were implemented using Python. The optimization of the hyperparameters were performed using the scikit-optimize package. The multipliers of the physics constraint terms of models 2 ($\mathbf{\rm GP^{\mathcal{L}_{\rm phy}}}$) and 4 ($\mathbf{\rm GP^{\rm upd,\mathcal{L}_{\rm phy}}}$) were chosen as $(\lambda^{\rm GP}_{\rm phy,1},\lambda^{\rm GP}_{\rm phy,2})=(50,50)$ based on a cross-validation test. The Automatic Relevance Determination (ARD) squared exponential function [43] was used as the covariance function for all the GP models.

The four DNN models (Models 5 to 8) were implemented using the Keras package [6] with Tensorflow in the backend. The hyperparameters of each model were tuned with grid search and the multipliers of the physics constraint terms in Model 6 ($\mathbf{\rm DNN^{\mathcal{L}_{\rm phy}}}$) and Model 8 ($\mathbf{\rm DNN^{\rm upd,\mathcal{L}_{\rm phy}}}$) were chosen as $(\lambda^{\rm DNN}_{\rm phy,1},\lambda^{\rm DNN}_{\rm phy,2})=(0.01,0.01)$ based on a cross-validation test. Fully-connected DNN models with 2 hidden layers and 5 neurons in each hidden layer were constructed. The Rectified Linear Unit (ReLU) activation function and Adam optimizer were used to perform stochastic gradient descent for 300 epochs in learning the model parameters. The dropout rate for the DNN models was chosen to be 0.05, for reasons as explained below.

The computational costs of different models for training and estimation of Sobol’ indices based on 5000 MC samples with a fixed number of experimental data ($n=39$) is given in Table. 1. Among the GP models, the time it takes for training as well as computation of Sobol’ indices using Models 2-4 is significantly greater than Model 1. The reason for the difference between the training time of $\mathbf{\rm GP}$ and $\mathbf{\rm GP^{\rm upd}}$ is the pre-training phase, where a large amount of physics input-output samples used. Whereas, the difference between the training time of $\mathbf{\rm GP}$ and $\mathbf{\rm GP^{\mathcal{L}_{\rm phy}}}$ is due to the inclusion of physics constraints, which makes it harder for the optimization to find optimal hyperparameters. Interestingly, the training time for the DNN models ranges from 20 to 55 sec on the same desktop computer as used for the GP model (Intel^® Xeon^® CPU E5-1660 v4$@$3.20GHz with 32 GB RAM and GPU NVIDIA Quadro K620 with 2 GB). Among the DNN models, the reasons for the increased training time of models 6-8 compared to Model 5 are the same as for the GP models. The Sobol’ index estimations based on 5000 samples take approximately 1-2 minutes for the DNN models (models 5-8), whereas the GP predictions take much longer because the covariance matrix needs to be stored and inverted.

In order to compare the accuracy of the eight different models, the models are trained with different amounts of experimental observations ($n=(5,10,15,20,30)$), and the remaining 9 observations are used to compute the errors; the root mean square error (RMSE) based on the 9 validation samples is reported as the accuracy measure for comparison. The mean and one standard deviation RMSE values for the four GP models and the four DNN models are shown in Tables 2 and 3 respectively, for different values of $n$. To further validate the accuracy of the models, the data set is divided into two subsets for cross-validation; $k$-fold cross-validation is performed by splitting the data set into sets for model training and cross-validation, and these sets are selected randomly $k=10$ different times. The average cross-validation accuracy of the models over the 10 folds (random shuffles) is assessed by evaluating the average RMSE (see Table 4).

The results in Tables 2, 3 and 4 show that the use of PIML strategies in DNN models improves the performance, and the improvement is relatively larger as the amount of observed data $n$ gets smaller. For $n$=20 and $n$=30, the basic DNN model is as accurate as the DNN models with the PIML strategies; whereas for smaller values of $n$, the DNN models with the PIML strategies are significantly more accurate. This clearly indicates the benefit of PIML over basic ML for the DNN models. However, the GP models incorporating PIML strategies do not show significant improvement. This is because the GP models have a much smaller number of parameters compared to the DNN models, and the achievable accuracy of any ML model is constrained by the number of model parameters (in addition to model form).

In addition to the increased number of parameters, the DNN models have additional advantages in terms of training epochs and dropout rate that affect the prediction accuracy. As mentioned earlier, the optimum number of training epochs was found to be 300. Regarding dropout rate, the RMSE values (based on 9 validation samples) of DNN models with MC dropout for different dropout rates are shown in Fig. 4. The lowest RMSE values for all four models are obtained between 0.005-0.05. On the other hand, the Sobol’ index estimates were found to be similar within this range of dropout rate, for each of the four DNN models. Therefore, a dropout rate of 0.05 was chosen for the reporting of GSA results, since a higher dropout rate results in smaller sub-networks, therefore a smaller number of parameters and faster training. By comparing Tables 2 and 3, it is seen that at the dropout rate of 0.05, the DNN models are more accurate compared to the GP models.

The Sobol’ indices estimates obtained using the two approaches of the second PIML strategy are compared in Figs. 5 and 6. The results show that both approaches converge to similar first-order and total effect sensitivity index estimates. Thus, in the rest of the paper the first approach of the second PIML strategy is used for both GP and DNN models, i.e., $\mathbf{\rm GP^{\rm upd}}$ and $\mathbf{\rm DNN^{\rm upd}}$.

The Sobol’ index computations with the GP models (1-4) are based on 5000 MC samples and 100 realizations of the Gaussian process. The effect of the number of experimental observations (used to train the GP models) on the first-order Sobol’ index estimates from the four GP models is illustrated in Fig. 7. And the total effect Sobol’ index estimates from the GP models for different numbers of experimental training data are shown in Fig. 8. The mean values of sensitivity estimates based on the GP model predictions are denoted with solid dots at a given number of observations $n$. The sensitivity results are reported for two input variables, nozzle temperature and nozzle speed; the output quantity of interest is porosity.

The 95% prediction intervals are represented with bars above and below the solid dots for the corresponding model. The bounds in the sensitivity estimates are calculated using the mean predictions based on the 100 realizations of the GP models. As expected, the prediction intervals decrease as increasing amounts of experimental data are used to train the GP models. Further, all four models converge to similar first-order and total effect sensitivity estimates for both printer nozzle temperature and speed. The relative individual contribution (at $n$=39) of nozzle speed to the variance of the porosity ($\approx$ 0.65) is greater than that of the nozzle temperature ($\approx$ 0.25) and their sum is $\approx$ 0.9. And the sum of their total effect indices (which capture parameter interactions) is slightly above 1.0, indicating that the interaction effect is small.

For further clarity regarding uncertainty, numerical values of the prediction bounds of first-order sensitivity estimates of nozzle temperature obtained using 100 realizations of the GP models are shown in Table 5. The results indicate that prediction intervals obtained using the PIML models 3 and 4 ($\mathbf{\rm GP^{\rm upd}}$ and $\mathbf{\rm GP^{\rm upd,\mathcal{L}_{\rm phy}}}$) converge to the bounds obtained using 39 number of observations for training the models faster than the first two models. Note that the upper 95% bound for Model 2 ($\mathbf{\rm GP^{\mathcal{L}_{\rm phy}}}$) converges to its final value (0.33) within 10 experimental observations.

The Sobol’ index computations with the DNN models (5-8) are based on 5000 MC samples and 100 stochastic forward passes through the networks. The calculated first-order sensitivity estimates of temperature and speed, $S_{T}$ and $S_{S}$ respectively, for different numbers of experimental observations (training data) $n=(5,10,20,30,39)$ are shown in Fig. 9. The distributions of sensitivity estimates are obtained using MC dropout predictions based on 100 stochastic forward passes through the networks for different number of observations. The mean values of sensitivity estimates are represented with solid dots and the 95% bounds are denoted with bars above and below the solid dots for the corresponding model. Similarly, the calculated total effect sensitivity estimates $S_{T}$ and $S_{S}$ for different values of $n$ are illustrated in Fig. 10. Similar to the GP results, the difference between the total effect and first-order indices of process inputs is negligible, which indicates that their interaction is not significant.

All DNN models converge to similar first-order and total effect sensitivity estimates for both inputs, and these values are consistent with the results obtained using GP models. For this problem with two inputs and one output, 39 experimental observations appear adequate to train even the basic DNN model to achieve similar performance as the other physics-informed models; in fact, the results are similar even at 20 observations. However, the superior accuracy of the physics-informed models becomes apparent if only a small number of experiments are available, say $n=5$ or $n=10$, as shown in Table 3 and discussed earlier in Section 4.1.4.

For further clarity regarding uncertainty, numerical values of the prediction bounds of first-order sensitivity estimates of nozzle temperature obtained using 100 forward passes through the DNN models are given in Table 6. The results show that prediction intervals obtained using all models show a similar trend. The 95% bounds for all models are significantly narrower than the ones obtained using the GP models. The uncertainty in the sensitivity estimates due to the DNN models is almost negligible when more than 10 number of observations are used to train the models.

The 95% prediction intervals (upper limit-lower limit) decrease as increasing amounts of experimental data are used to train the DNN models. However, the prediction intervals obtained using the DNN models are much smaller than the ones obtained using the GP models since the DNN models has more degrees of freedom that can be optimized. For example, for $n=39$, the prediction interval width is 0.02-0.03 for the DNN models, whereas it is 0.12-0.15 for the GP models. In addition to the larger number of parameters, the number of training epochs, which is the number of complete passes through a batch of training dataset, is optimized for the DNN models, thus maximizing their prediction accuracy. The prediction accuracy is further improved by choosing the appropriate dropout rate as discussed earlier. The number of parameters to be learned reduces with the use of dropout, which helps with regularization and prevents ill-conditioning. Further, the number of training epochs is also observed to affect both the accuracy and uncertainty of the sensitivity estimates. It was found that at a smaller number of epochs, the model was not fully trained resulting in underfitting, leading to larger bias and variance in the prediction. As the number of epochs was increased, both the bias and the variance were reduced.

In order to further illustrate the effectiveness of the proposed method and how different combinations of methods affect the accuracy and computational effort, a more complex problem involving eleven random variables is considered. The data set used by Karpatne et al. [26] for lake temperature modeling is considered, at Lake Mendota in Wisconsin, USA.

The overall data for Lake Mendota comprised of 13,543 temperature observations from 30 April 1980 to 02 Nov 2015. The modeling uses used 11 meteorological drivers as input variables. The original physics-based model, referred to as the General Lake Model (GLM), modeled the lake temperature by performing 1-D analysis (along depth) considering a variety of lake variables, and produced a total of 662,781 training samples (input-output). The physical non-linear relationship between temperature, depth and density of water is used to evaluate the physical violations across every consecutive depth-pair and time step. A more detailed description of the physics-based loss function and data can be found in the work by Karpatne et al. [26].

The computational costs of different models for training and estimation of Sobol’ indices based on 5000 MC samples with a fixed number of experimental data ($n=1000$) is given in Table. 7. The desktop computer used for the illustrative example 1 (Intel^® Xeon^® CPU E5-1660 v4$@$3.20GHz with 32 GB RAM and GPU NVIDIA Quadro K620 with 2 GB) is used for the computations of Sobol’ indices. Among the GP models, the time it takes for training as well as computation of Sobol’ indices using Models 2-4 is similar to Model 1. The Sobol’ index estimations based on 5000 samples take approximately 1-3 minutes for the DNN models (Models 5-8). The training time and the computation of Sobol’ indices using DNN models are significantly less than Models 1-4, which use the GP models, since the GP predictions require the covariance matrix to be inverted.

The effect of the number of experimental observations on the first-order and total effect Sobol’ index estimates from the four GP models of the lake temperature modeling is illustrated in Figs. 11 and 12.

The 95% prediction intervals are represented with bars above and below the solid dots for the corresponding model. The prediction intervals decrease as increasing amounts of experimental data are used to train the GP models. In addition, all four models converge to similar first-order and total effect sensitivity estimates for all input variables.

The Sobol’ index computations with the DNN models (5-8) are based on 5000 MC samples and 100 stochastic forward passes through the networks. The calculated first-order sensitivity estimates of eleven input variables, for different numbers of experimental observations (training data) $n=(100,250,500,750,1000)$ are shown in Figs. 13 and 14.

All DNN models converge to similar first-order and total effect sensitivity estimates for all eleven inputs, and these values are consistent with the results obtained using GP models. The total effect sensitivity estimates obtained using the basic GP and DNN models for most of the input variables converge to a slightly different value than the ones obtained using the proposed PIML models. Note that the mean estimates of the PIML models converge more smoothly than the basic ML models.

The results show that prediction intervals obtained using all models show a similar trend. The 95% bounds for all DNN models are significantly narrower than the ones obtained using the GP models. The results indicate that the mean estimates and prediction intervals obtained using the PIML models 4 and 8 ($\mathbf{\rm GP^{\rm upd,\mathcal{L}_{\rm phy}}}$ and $\mathbf{\rm DNN^{\rm upd,\mathcal{L}_{\rm phy}}}$) converge to the mean estimates and bounds obtained using 1000 number of observations for training the models faster than the other models. Note that the mean estimates for Models 4 and 8 ($\mathbf{\rm GP^{\rm upd,\mathcal{L}_{\rm phy}}}$ and $\mathbf{\rm DNN^{\rm upd,\mathcal{L}_{\rm phy}}}$) converge to their final values within 250 experimental observations for almost of the variables. The prediction bounds of DNN models do not change significantly as more experimental observations are used to train the models. However, the prediction bounds get smaller as more number of training data is used to train the GP models. Moreover, the prediction bounds of physics-informed models are slightly narrower than the basic GP model. For further clarity regarding uncertainty, numerical values of the prediction bounds of first-order sensitivity estimates of $X_{1}$ (day of year) obtained using 100 realizations of the GP models and 100 forward passes through the DNN models are shown in Table 8.

The results of these numerical examples could be summarized as follows:

• •

  The GP models required more computational effort than the DNN models, both in training and prediction.

• •

  The DNN models were able to achieve higher accuracy and lower uncertainty in prediction, due to the optimization of dropout rate and number of training epochs.

• •

  The physics-informed ML models are able to achieve higher prediction accuracy than the basic ML models, especially when the amount of available experimental data is small.

Overall, for these numerical examples, the DNN models gave higher accuracy and lower uncertainty in the GSA results than the GP models, and also required less computational effort both in training and prediction. However, the first numerical example consisted of only two inputs and a single output. As the number of inputs and outputs increase with the second numerical example, all the ML models considered above were challenged w.r.t. adequacy of training data, computational effort in training and prediction, and the accuracy and uncertainty of the sensitivity estimates. The results showed that the mean estimates of the PIML Models 4 and 8 ($\mathbf{\rm GP^{\rm upd,\mathcal{L}_{\rm phy}}}$ and $\mathbf{\rm DNN^{\rm upd,\mathcal{L}_{\rm phy}}}$) for both numerical examples converge to the mean estimates obtained using the largest number of observations (i.e., 39 and 1000 observations for the first and second numerical example, respectively) for training the models faster than the other models.