Model uncertainty is a significant source of epistemic uncertainty that affects the prediction of a multidisciplinary system. In order to achieve a reliable design, it is critical to ensure that the disciplinary/subsystem simulation models are trustworthy, so that the aggregated uncertainty of system quantities of interest (QOIs) is acceptable. Reduction of model uncertainty can be achieved by gathering additional experiments and simulations data; however, resource allocation for multidisciplinary design optimization (MDO) and analysis remains a challenging task due to the complex structure of the system, which involves decision makings about where (sampling locations), what (disciplinary responses), and which type (simulations versus experiments) for allocating more resources. Instead of trying to concurrently make the above decisions, which would be generally intractable, we develop a novel approach in this paper to break the decision making into a sequential procedure. First, a multidisciplinary uncertainty analysis (MUA) is developed to identify the input settings with unacceptable amounts of uncertainty with respect to the system QOIs. Next, a multidisciplinary statistical sensitivity analysis (MSSA) is developed to investigate the relative contributions of (functional) disciplinary responses to the uncertainty of system QOIs. The input settings and critical responses to allocate resources are selected based on the results from MUA and MSSA, with the aid of a new correlation analysis derived from spatial-random-process (SRP) modeling concepts, ensuring the sparsity of the selected inputs. Finally, an enhanced preposterior analysis predicts the effectiveness of allocating experimental and/or computational resource to answer the question about which type of resource to allocate. The proposed method is applied to a benchmark electronic packaging problem to demonstrate how epistemic model uncertainty is gradually reduced via resource allocation for data gathering.

## Introduction

With the increasing accuracy and complexity of simulation models, the cost of running high-fidelity computer simulations becomes more and more prohibitive in simulation-based design. Accordingly, emulators (also known as metamodels or surrogate models) have been widely used as a replacement of original simulation models and are especially useful in searching for the optimal design. However, due to the lack of simulation/experimental data, the accuracy of emulators can be very poor, which introduces *epistemic* model uncertainty that, in contrast to the *aleatory* uncertainty from natural/physical randomness, can be reduced by collecting additional data. We refer to the process of planning additional experiments/simulations for reducing epistemic model uncertainty as a *resource allocation* problem. Resource allocation is even more important in MDO, which, unlike traditional design optimization problems, is notoriously complicated due to the fact that it requires analyses in multiple disciplines and/or involves a number of subsystems and components. In particular, resource allocation for MDO requires an optimal scheme to distribute resources to different disciplines and to enhance the prediction of system QOIs. However, this is not an easy problem, because the coupling relationship between different disciplines leads to a nested structure and substantial complexity.

Model uncertainty, as a major source of epistemic uncertainty in multidisciplinary systems, stems from either having a limited knowledge of the true reality and/or using a simplified model that does not satisfactorily represents the reality (represented as “model discrepancy” [1]), or from a lack of computer simulation data (termed “code uncertainty” [1] or “interpolation uncertainty” [2]) when replacing the expensive computer model using a surrogate model (metamodel). Even though the sources of these two types of uncertainty are different, they all belong to epistemic model uncertainty due to the lack of data—model discrepancy uncertainty due to the lack of experimental data and interpolation uncertainty due to the lack of computer simulation data. Both sources of epistemic uncertainty can be addressed simultaneously using the Bayesian Gaussian process approach presented by Kennedy and O'Hagan [1]. The epistemic model uncertainty of system QOIs is an aggregation of the uncertainties in the disciplinary subsystems, and the purpose of resource allocation is to reduce the aggregated prediction uncertainty of system QOIs. Nevertheless, it is not straightforward to provide a good resource allocation scheme, due to the amount of decisions involved and their dynamic nature. The decisions to make include:

- (1)
*Where*in the input space of a multidisciplinary system should we allocate additional resources? (That is, what are the input locations where system QOIs are most influenced by epistemic model uncertainty?) - (2)
To

*what*disciplinary response(s) shall we allocate more resources? - (3)
*Which*type of resources shall we allocate? Is adding more simulations sufficient, or should we conduct more physical experiments?

The first question requires uncertainty quantification of disciplinary subsystems and an efficient way to propagate their impact to the system QOIs. Although there has been extensive research on uncertainty propagation of a single discipline [1,3–8], or a multidisciplinary system but with only aleatory uncertainty (e.g., using Monte Carlo simulation [9–11], Taylor series expansion [12–15], decomposition strategies [16–22], etc.), existing work for multidisciplinary systems with epistemic model uncertainty is scarce. The problem is often oversimplified by either predefining model bias [23–27] or assuming a simple form of interdisciplinary couplings (restricted to feed-forward couplings [28], for example). Furthermore, even after developing a convenient uncertainty propagation scheme, one must still decide where to choose the input sample locations so that they sufficiently cover all regions with large epistemic model uncertainty. Consequently, the second question requires an approach to distinguish the relative contributions (which can vary over the input space) of the different disciplinary simulation models on the uncertainty of system QOIs. Statistical sensitivity analysis (SSA) is commonly used in engineering design to gain knowledge of complex model behavior and to facilitate designers' decision making on where to spend engineering effort [29]. However, SSA of model uncertainty, which is different from SSA of random variables in the more traditional literature, has rarely been studied [30–32]. The difficulty lies in the fact that uncertain models are essentially functional inputs rather than scalars, and hence, their SSA would be fundamentally different than traditional routines. Allaire et al. [30] developed a pure sampling-based approach, whose applicability is practically limited to single-disciplinary systems and hierarchical systems without strong couplings due to its computational cost. SSA for a fully coupled multidisciplinary system was never studied by other researchers. The third question is perhaps the least studied of the three. Previous works focused on allocating a single type of resource, typically experimental (e.g., Ref. [28]). Indeed, relative to conducting simulations, conducting physical experiments may be a more direct and effective way to learn physical reality, but it is also typically more expensive. For that reason, there exists a tradeoff between most reducing uncertainty (by conducting experiments) versus least increasing computational budget (by conducting simulations).

On top of these three individually challenging questions, the complexity of concurrently answering them adds another layer of difficulty. To simplify the problem, the pioneering work by Sankararaman et al. [28] had to assume no coupling relationship between disciplines, and that the model uncertainty is purely due to some unknown model parameters, which makes SSA much easier.

To manage the complexity in the decision-making process, this paper presents a new strategy and mathematical framework for resource allocation in a multidisciplinary system, one that breaks the decision making into a sequential procedure. An MUA method [33,34] we developed earlier is applied to make the first decision (i.e., *where*). A multidisciplinary sensitivity analysis (MSSA) [31,32], which we developed to handle uncertain *functional responses* as the stochastic inputs in the sensitivity analysis (the stochastic inputs are scalar variables in the traditional sensitivity analysis literature), is applied to investigate the relative contributions of disciplinary responses to the aggregated uncertainty of system QOIs for the second decision (i.e., *what*). A type of correlation analysis derived from SRP modeling concepts is also developed as an aid to ensure the sparsity of the selected input settings. Finally, an enhanced preposterior analysis predicts the effectiveness of allocating experimental or computational resource, or a combination of both, to make the last decision (i.e., *which*). The remainder of the paper is organized as follows. Section 2 reviews the SRP-based model bias correction approach to quantify epistemic model uncertainty and its extension in multidisciplinary systems. Section 3 details the proposed framework for resource allocation, which is applied in Sec. 4 to a benchmark electronic packaging problem to illustrate how the uncertainty of system QOIs is reduced over the input space via a few iterations of resource allocation. Concluding remarks are provided in Sec. 5.

For clarity of the discussions in this paper, we make the following remarks:

- (1)
The objective of this paper is to improve the

*global*modeling capability of a multidisciplinary system, such that the epistemic model uncertainty of system QOIs is acceptable over the input space. A globally accurate model can be useful in many design activities other than direct optimization; for example, the simulation/emulation results of a model can be passed into a classification analysis for designers to understand the potentially feasible/infeasible input settings and to be more focused on the “regions of interest.” This line of research [35,36] is different from the objective-oriented sampling [37–39] which specifically targets a*local*region where potential optimal design solution lies. - (2)
When referring to

*resources*, we mean*experiments*and*simulations*. We use “experiments” in a broad context: it can be real physical experiments, or extremely high-fidelity simulations that (almost) can represent the reality. “Simulations” in this paper generally refer to results from physics-based models that are less accurate than experiments. - (3)
*Emulators*are nonphysics-based models that are built upon data from experiments and simulations. They are much less costly than direct simulations. In this work, an SRP modeling technique is used to create disciplinary emulators. - (4)
Experiments and simulations are conducted for individual disciplinary responses, but not for the system QOIs (which we consider practically unaffordable). Any full systemwise analysis is conducted using emulators.

- (5)
The framework proposed in this paper is targeted at a multidisciplinary problem, although it can be simplified to work for a single discipline as well.

## Review of Model Bias Correction in Multidisciplinary Systems

Model uncertainty is a significant type of epistemic uncertainty in a multidisciplinary system. To quantify model uncertainty, data from the computer simulations are first collected and then compared with the experimental data. Based on the differences in the response values over the two datasets, the original simulation model is “updated” to a more accurate emulator; meanwhile, the model uncertainty is quantified. We refer to this process as *model bias correction*. Subsequently, an extra experimental dataset is provided to validate the emulator. Model bias correction is conducted iteratively by refining the model and/or collecting additional data, until the emulator is trustworthy. In this section, we briefly review the SRP-based model bias correction technique and its application in multidisciplinary systems.

where the superscript “*e*” denotes the experimental response of *y*(**x**), and the superscript “*m*” denotes the simulation model of the same response. *δ*(**x**) is the discrepancy function that represents the model bias, and *ε* ∼ $N$(0,*λ*) is a random error accounting for the experimental variability, assumed to be normally distributed with unknown variance *λ*. The situation where *ε* does not have zero mean (i.e., a so-called experimental bias or “systematic error”) would be more complicated and is beyond the scope of this paper.

*M*simulation response observations ${xim,ym(xim):i=1,\u2026,M}$ and a set of

*N*experimental response observations ${xie,ye(xie):i=1,\u2026,N}$ have been collected. Since $xim(i=1,\u2026,M)$ and $xie(i=1,\u2026,N)$ are not necessarily the same input settings, and likewise, $ym(xim)$ and $ye(xie)$ may not be directly compared to obtain the discrepancy function, an SRP modeling technique is adopted to bridge the gap and enable the comparison. An SRP can be viewed as a collection of random variables distributed over the input space that collectively have certain multivariate spatial correlation structure. Specifically, Gaussian process (GP) modeling, which assumes the random variables to follow a multivariate normal distribution, has emerged as one of the most popular metamodeling approaches due to its ability to capture functional nonlinearity, perfectly interpolate deterministic surfaces, etc. [40,41]. In this study, we adopt GP modeling and model the simulation and the discrepancy function as GPs such that

In the preceding, **h*** ^{m}*(

**x**) and

**h**

*(*

^{δ}**x**) denote vector-valued functions whose elements are some user-specified regression basis functions (e.g., constant, linear, quadratic, etc.), and

**β**

*and*

^{m}**β**

*are two vectors of regression coefficients associated with*

^{δ}**h**

*(*

^{m}**x**) and

**h**

*(*

^{δ}**x**), respectively. Their products,

**h**

*(*

^{m}**x**)

^{T}**β**

*and*

^{m}**h**

*(*

^{δ}**x**)

^{T}**β**

*, construct the prior mean functions of*

^{δ}*y*(

^{m}**x**) and

*δ*(

**x**).

*σ*and

_{m}*σ*are the prior standard deviations, and

_{δ}**ω**

^{m}^{ }= $[\omega 1m,\omega 2m,\u2026,\omega pm]T$ and

**ω**

^{δ}^{ }= $[\omega 1\delta ,\omega 2\delta ,...,\omega p\delta ]T$ are the spatial correlation parameters used to capture the nonlinearity of the functions.

*V*(

^{m}**x**,

**x**′) and

*V*(

^{δ}**x**,

**x**′), which depend on

*σ*,

_{m}*σ*,

_{δ}**ω**

*, and*

^{m}**ω**

*, are the covariance functions of*

^{δ}*y*(

^{m}**x**) and

*δ*(

**x**), respectively.

where $d={ym(x1m),\u2026,ym(xMm),ye(x1e),\u2026,ye(xNe)}T$ is the collected response data from simulations and experiments, and **β** = [(**β*** ^{m}*)

*, (*

^{T}**β**

*)*

^{δ}*]*

^{T}*is the collection of regression coefficients.*

^{T}**H**is a matrix generated by

**h**

*(*

^{m}**·**) and

**h**

*(*

^{δ}**·**), and

**V**is generated by

_{d}*V*(·,·) and

^{m}*V*(·,·) (see Refs. [8] and [41] for full derivation and discussion).

^{δ}**ϕ**= {

**β**

*,*

^{m}*σ*,

_{m}**ω**

*,*

^{m}**β**

*,*

^{δ}*σ*,

_{δ}**ω**

*, λ} denote the collection of unknown parameters of the GP models, referred to as*

^{δ}*hyperparameters*. The key to estimating the discrepancy function and subsequently enabling a better model prediction is to estimate these hyperparameters. They are usually estimated via

*maximum likelihood estimation*, i.e., maximizing the likelihood function

*p*(

**ϕ**|

**d**). Subsequently, we can obtain a prediction of

*y*(

^{e}**x**) at any untested design of interest

**x**. The relationship between

*y*(

^{e}**x**) and its prediction can be written as

*δ*(

**x**).

*Z*(

**x**) is a zero-mean GP, i.e.,

*Z*(

**x**)|

**x**∼ $N$(0, $\sigma Z2(x)$), which accounts for the uncertainty in the prediction. $\sigma Z2(x)$ is also known as the mean square error of the mean prediction $y\u0302e(x)$. The merit of using GP models is that we are able to obtain analytical equations for calculating $y\u0302e(x)$ and $\sigma Z2(x)$. For example,

Figure 1 illustrates a typical result from GP-based model bias correction, in which there is originally a significant model bias when comparing the simulation data (the triangles) and the experimental data (the circles), especially on the left and right boundaries of the *x* domain. Using GP modeling, we are able to obtain an updated model (the solid purple curve) that matches the experimental data, while maintaining the general functional trend learned from the simulation data. Besides, GP modeling enables the quantification of interpolation uncertainty at the locations that have not yet been simulated, represented by a prediction interval (PI) (the shaded region) calculated from $\sigma Z2(x)$.

GP-based model bias correction can be applied to multidisciplinary systems. A notional multidisciplinary system is depicted in Fig. 2. Suppose it involves a total of ND disciplines associated with a collection of individual disciplinary input variables **x**_{ind} = {**x*** _{i}*:

*i*= 1,…, ND} and some input variables

**x**

*that are shared across at least two disciplines. The disciplines are coupled via linking variables*

_{s}**u**

*(*

_{ij}*i*,

*j*= 1,…, ND), which serve as output from the

*i*th discipline and input to the

*j*th discipline. If both

**u**

*and*

_{ij}**u**

*are nonempty, the relation between the*

_{ji}*i*th and

*j*th disciplines is referred to as

*feedback*coupling, and otherwise,

*feed-forward*coupling. We denote the collection of all linking variables that are output from and input to the

*i*th discipline as

**u**

*. and*

_{i}**u**.

*, respectively, i.e.,*

_{i}**u**

*. = {*

_{i}**u**

*:*

_{ij}*j*= 1, …, ND

*, j*≠

*i*} and

**u**.

*= {*

_{i}**u**

*:*

_{ji}*j*= 1, …, ND

*, j*≠

*i*}. The system QOIs

**y**

_{sys}depend on individual disciplinary responses

**y**

_{ind}= {

**y**

*:*

_{i}*i*= 1, …, ND} collected from the responses of ND disciplines, through the system analysis model.

*i*th discipline, the simulated response

**u**

*. may have uncertainty, which will introduce uncertainty to*

_{i}**u**.

*and further influence*

_{i}**u**

*. itself. Such “convolution” accumulates uncertainty in the disciplinary level and passes it on to the disciplinary outputs*

_{i}**y**

*. Another complication is that the simulation model for*

_{i}**y**

*may also involve epistemic model uncertainty. That is, even if*

_{i}**u**.

*has no uncertainty, simulating the response*

_{i}**y**

*still could involve uncertainty. Based on Eq. (4), we present the predictions of all the linking variables and the disciplinary output variables, respectively, as (for*

_{i}*i*= 1, 2,…, ND)

Assessing the cumulative uncertainty of system QOIs calls for an uncertainty propagation approach. Our recent development of a GP-based MUA method [33,34] addresses the issues of efficient uncertainty propagation in multidisciplinary systems, which is generically suitable for situations with both aleatory uncertainty (from input variables **x**_{ind} and **x*** _{s}*) and epistemic uncertainty (from models). A version of the method that only considers epistemic uncertainty (which is the main focus of this paper) will be provided in Sec. 3.1.

## Resource Allocation for Uncertainty Reduction in Multidisciplinary Systems

In this section, we present in detail the proposed resource allocation approach for reducing epistemic model uncertainty in a multidisciplinary system. An overview is provided in Fig. 3. In the preliminary step, model bias correction is conducted based on some existing data to construct emulators for disciplinary responses (linking the input variables **u*** _{ij}* and

**x**

*to the disciplinary outputs*

_{j}**y**

*, for each*

_{i}*i*). In step 1, a space-filling strategy is used to explore the input space {

**x**

_{ind},

**x**

*} by generating sufficient sample points and identifying the locations with unacceptable amounts of uncertainty with respect to the system QOIs*

_{s}**y**

_{sys}. The MUA method is applied to assess the aggregated epistemic model uncertainty of

**y**

_{sys}at each sample point (Sec. 3.1). An MSSA then separates the impact of epistemic model uncertainty from different responses at each sample point (Sec. 3.2). Decisions about resource allocation (

*where*,

*what*, and

*which*) are made in sequence in step 2, where we integrate the information from MUA, MSSA, and a preposterior analysis. Finally in step 3 as resources are allocated, another round of model bias correction is performed.

**y**

_{sys}is an indicator of the amount of uncertainty, we specify the following stopping criterion in this paper:

The numerator is the standard deviation of **y**_{sys} (at any given input setting {**x**_{ind}, **x*** _{s}*}), and the denominator is the mean of the absolute value of

**y**

_{sys}over the input space. Thus,

*γ*is an intuitive measure of the prediction uncertainty that resembles the “coefficient of variation” in traditional statistics. The user can choose other variance-based stopping criteria, independent from the rest of procedure presented.

### Exploring the Input Space and Evaluating the Impact of Epistemic Model Uncertainty.

In step 1, to find where in the input space is most influenced by epistemic model uncertainty, a straightforward treatment is to apply a space-filling strategy (for example, Monte Carlo uniform sampling, Latin hypercube sampling (LHS), Halton sequence, etc.) and generate *N _{s}* samples over the input space: $(xind(k),xs(k))=(x1(k),\u2026,xND(k),xs(k)),k=1,\u2026,Ns$. Next, for each sample point, the GP-based MUA method [33,34] is applied to evaluate the mean and variance of system QOIs

**y**

_{sys}. The merit of the MUA method is that it utilizes the structure of GP emulators and enables an efficient analytical uncertainty analysis. The following discussion in this subsection is for every sample input that is generated, and hence for notational simplicity, we omit the superscript “(

*k*).”

**x**

_{ind},

**x**

*), a first-order Taylor expansion of $u\u0302i\xb7e(\xb7,\xb7,\xb7)$ and $y\u0302ie(\xb7,\xb7,\xb7)$ in Eqs. (7) and (8) is applied, which leads to (for*

_{s}*i*= 1,…, ND)

**μ**

*.,*

_{ui}**μ**

*.*

_{u}*, and*

_{i}**μ**

*denote the means of*

_{yi}**u**

*.,*

_{i}**u**.

*, and*

_{i}**y**

*, respectively. The derivatives on the right-hand side, while not explicitly written out, are taken at the mean of the variables, i.e., (*

_{i}**x**

_{ind},

**x**

*,*

_{s}**μ**

*). Therefore, the means and covariance matrices of the linking variables (i.e., all*

_{u}**u**

*'s) and the disciplinary outputs (i.e., all*

_{ij}**y**

*'s) can be estimated by (for*

_{i}*i*= 1,…, ND)

**Σ**

*,*

_{u}**Σ**

*,*

_{Zu}**Σ**

*, and*

_{y}**Σ**

*denote the covariance matrices of*

_{Zy}**u**

*,*

^{e}**Z**

*,*

_{u}**y**

*, and*

^{e}**Z**

*, respectively. The matrices*

_{y}**A**and

**B**are

Since GP-based model bias correction provides analytical formulas to calculate $u\u0302i\xb7e(\xb7,\xb7,\xb7)$ and $y\u0302ie(\xb7,\xb7,\xb7)$ (similar to Eq. (5)), the entries of matrices **A** and **B** also have analytical forms, which enables an efficient calculation.

Detailed derivations can be found in Refs. [33] and [34]. After deriving the means and covariance matrices of the disciplinary responses **y**_{ind}, the derivation of the mean and variance of the system QOIs **y**_{sys} can be treated as a single-disciplinary uncertainty propagation problem and solved using any conventional method. From the variance of **y**_{sys}, it can be easily discovered which regions have large uncertainty.

### Separating the Impact of Epistemic Model Uncertainty From Disciplinary Responses.

Apart from assessing the aggregated impact of epistemic model uncertainty on **y**_{sys} over the *N _{s}* sample points, one can actually separate the total impact into the contributions from different responses (

**u**

*'s and*

_{ij}**y**

*'s) using SSA. This information is particularly important for resource allocation in a multidisciplinary system, as it helps decision makers understand what the most important factors are that influence the system performance. Moreover, it will be shown in Sec. 3.3 that it also helps in choosing the most appropriate input settings that are most in need of additional resources.*

_{i}The most widely used measure of variability in SSA is variance, because it is usually interpretable for practitioners. Building on the Sobol’ method [42,43], which is a variance-based Monte Carlo method for aleatory input variability, we developed a more comprehensive MSSA in Refs. [31] and [32] that will assess the impact of aleatory and epistemic uncertainties. Different from the Sobol’ method, MSSA has the capability of measuring the impact of *functional responses* that serve as stochastic inputs in the sensitivity indices (the stochastic inputs are scalar variables in the traditional Sobol’ method, although the output is a functional response). Particularly, we use *local* sensitivity analysis in this subsection, which focuses on examining the impact of uncertainties from different sources at a single input setting.

*Y*can be decomposed into the contributions of a set of random input variables

*X*

_{1},…,

*X*. The

_{p}*main sensitivity index*(MSI) and the

*total sensitivity index*(TSI) of a particular input

*X*, which measure the main effect and the overall effect of

_{i}*X*on the output

_{i}*Y*, respectively, are given by

**X**

_{∼}

*denotes all the input variables except*

_{i}*X*. The subscript in “Var” and “ $E$ ” denotes the source of uncertainty; for example, $VarXi(\u22c5)$ is the variance of the variable in the brackets due to the variation of

_{i}*X*. Both MSI and TSI are between 0 and 1, with 0 indicating no effect and 1 the most dominant effect. The same idea can be extended to examining the relative contributions of epistemic model uncertainty. If we view any

_{i}*Z*∈{

_{l}**Z**

*.,*

_{ui}**Z**

*:*

_{yi}*i*= 1,…, ND} (the random term that quantifies epistemic model uncertainty of a particular response in Eqs. (7) and (8)) in a similar way as

*X*is treated in Eqs. (16) and (17), then the MSI and TSI of

_{i}*Z*on the system QOIs

_{l}**y**

_{sys}can be defined as

where **Z**_{∼}* _{l}* = {

**Z**

*.,*

_{ui}**Z**

*:*

_{yi}*i*= 1,…, ND}\

*Z*. The Sobol’ variance decomposition assumes that the uncertainty sources are uncorrelated, which is reasonable in our case because the epistemic uncertainty of the individual disciplinary models (i.e.,

_{l}*Z*'s) is quantified separately and not correlated with each other.

_{l}It is worth noting that although Eqs. (18) and (19) resemble the standard Sobol’ indices (Eqs. (16) and (17)), the calculation is much more complex due to the fact that the **Z**'s are the stochastic “input variables” in our indices, but they are really stochastic functional responses over the input variables (i.e., random processes rather than random variables). In standard Sobol’ indices, the output is a functional response over the stochastic inputs, but the stochastic inputs are scalar variables and not functional responses themselves. For example, a sampling-based treatment to calculate the MSI of *Z _{l}* as defined in Eq. (18) is a trilevel procedure: in the outer level, one needs to generate multiple Monte Carlo realizations of

*Z*(which are, by the nature, realizations of a function); in the intermediate level, for each realization of

_{l}*Z*, one needs to generate realizations of

_{l}**Z**

_{∼}

*(realizations of a set of functions); and finally in the inner level, given the realizations of all*

_{l}**Z**'s, computations are needed to evaluate $EZ\u223cl(ysys|Zl)$. Therefore, the overall computation could be extremely expensive. Likewise for the process of calculating the TSI of

*Z*.

_{l}The MUA method developed in Sec. 3.1 significantly speeds up the analysis. By applying Eqs. (12) and (13), $EZ\u223cl(ysys|Zl)$ is directly obtained without having to sample **Z**_{∼}* _{l}*, which eliminates the need of the outer level aforementioned. If necessary, one can further simplify the analysis and obtain rough approximations of MSI and TSI by linearizing $EZ\u223cl(ysys|Zl)$ and $EZl(ysys|Z\u223cl)$ as $[ysys|Z\u223cl\u22610,Zl]$ and $[ysys|Zl\u22610,Z\u223cl]$, respectively, i.e., fixing the randomly varying

**Z**terms to their zero means, which then eliminates the need of the intermediate level as well. Such linearization would certainly lead to a less accurate sensitivity assessment, but it is very useful for a fast analysis in a large complex system.

### Choosing Input Settings and Corresponding Responses for Resource Allocation.

To decide where to allocate more resources, a simple approach is to choose a few samples (up to the budget) with the largest variance of **y**_{sys} among the *N _{s}* samples generated in Sec. 3.1. However, this is not the best approach, since multiple sample points may happen to fall inside the same region that has a large epistemic model uncertainty, especially when

*N*is large. In such case, allocating resource to a single point inside this region should be sufficient for reducing the uncertainty. A good decision-making criterion is needed to detect whether two points with large uncertainty are “sufficiently far away” from each other, so that our final selection of samples would be widely spread over the input space without dense clustering. In this subsection, we propose a form of correlation analysis that utilizes information from the GP-based correlation, which measures the spatial proximity between two input settings, to select the best subset of input settings to allocate more resources.

_{s}*N*samples we generated in Sec. 3.1. We start by calculating their sensitivity indices and detecting the main contributing sources of uncertainty, respectively, at these two input locations. If their major sources of uncertainty are different, i.e., if $Var[ysys(xind(k),xs(k))](k=1,2)$ are mostly induced by model uncertainty from two different responses, then it is concluded that $(xind(1),xs(1))$ and $(xind(2),xs(2))$ are “far away” from each other and reside in different regions. On the other hand, if $Var[ysys(xind(k),xs(k))](k=1,2)$ turn out to be induced by uncertainty of the same response (for example, some response

_{s}*L*(

**x**

*,*

_{L}**x**

*,*

_{s}**u**.

*)∈{*

_{L}**u**

*.,*

_{i}**y**

*:*

_{i}*i*= 1,…, ND} associated with model uncertainty term

*Z*), then the two points may or may not be far away. We need to proceed to evaluate their correlation coefficient in the GP emulator of

_{l}*L*, defined as

In the denominator, $\sigma mL2$, $\sigma \delta L2$, and $\lambda L$ are the prior variances of the simulation model, the discrepancy function, and the experimental variance of response *L*, respectively; their sum, according to Eq. (2), is the prior variance of experimental response *L ^{e}* at $(xind(k),xs(k))(k=1,2)$. In the numerator,

*V*(·,·) and

^{mL}*V*(·,·) are the covariance functions of the simulation model and the discrepancy function of

^{δL}*L*, respectively, which can be directly calculated from Eq. (2); their sum is the prior covariance of the two sample points. $\mu u\ufffdL(k)$$(k=1,2)$ are the estimated values of linking variables (as inputs to the response

*L*). Therefore,

*ρ*

_{1,2}is the GP-based correlation coefficient, which is an indicator of how much impact the two input settings have over each other; in other words, how likely it is that they would yield a similar value of response and/or share a similar amount of uncertainty. It can be seen that this type of correlation coefficient is easy to compute and yet comprehensively combines the information of the absolute spatial distance between two sample points, the roughness of the response (quantified by

**ω**

*and*

^{m}**ω**

*in Eq. (2)), and the general trend of the response (quantified by*

^{δ}**β**

*,*

^{m}**β**

*,*

^{δ}*σ*, and

_{m}*σ*in Eq. (2)). If

_{δ}*ρ*

_{1,2}is close to 1 (or larger than a user-defined limit

*β*), the two sample points are considered to be in the same local region, and we conclude that allocating resource to one of them should be sufficient.

The Appendix gives a pseudocode for choosing the proper sample points for resource allocation. Furthermore, the developed method allows us to decide at each chosen sample point what response to allocate more resource, by picking the response *L* that mostly influences **y**_{sys} at this point.

### Deciding the Type of Resource to Allocate for Each Chosen Sample Point.

Finally, the one remaining task is to decide which type of resource, simulation or experiment, to allocate for each chosen input setting and the corresponding chosen response at this input setting. To address this, in this subsection, we propose an enhanced *preposterior analysis* approach (see Fig. 4, based upon a similar idea in the earlier work [44,45]), which, prior to allocating the resources, can *predict the posterior variance* of **y**_{sys} if the resources are actually gathered. Different from the existing work in the literature that focuses only on a single type of resource (usually experimental), the preposterior analysis proposed in this section provides insight for any combination of simulations and experiments.

Following the analyses in Sec. 3.3, suppose we have already made the decision to allocate more resources to response *L* at input locations $(xind(k),xs(k),\mu u\u22c5L(k))(k=1,\u2026,NL)$, and to response *L*′ at input locations $(xind(k),xs(k),\mu u\u22c5L\u2032(k))(k=NL+1,\u2026,NL+NL\u2032)$, etc. The preposterior analysis first generates hypothetical simulation and/or experimental data and then evaluates the influence on the posterior variance of the system QOIs. The detailed procedure is as follows:

- (1)
Suggest a resource allocation plan and generate hypothetical data.

To assess the impact of conducting experiments at $(xind(k),xs(k),\mu u\u22c5L(k))(k=1,\u2026,NLe)$ and conducting simulations at $(xind(k),xs(k),\mu u\u22c5L(k))(k=NLe+1,\u2026,NL)$ for response

*L*, hypothetical experimental and simulation data are generated accordingly at these locations. A preferred way is to generate them simultaneously considering their correlation. It should be noted that*L*(i.e., the experimental response of^{e}*L*) at $(xind(k),xs(k),\mu u\u22c5L(k))(k=1,\u2026,NLe)$ and*L*(i.e., the simulation response of^{m}*L*) at $(xind(k),xs(k),\mu u\u22c5L(k))(k=NLe+1,\u2026,NL)$ follow a jointly normal distribution that can be determined from the already estimated hyperparameters of the GP models. The equation for the jointly normal distribution is similar to Eq. (3), i.e.,**d**∼ $N$(**Hβ**,**V**), where_{d}**d**here should be interpreted as the hypothetical data to be generated, and**V**is a covariance matrix whose entries are the covariances between selected input settings and/or between_{d}*L*and^{e}*L*. Therefore, a set of hypothetical experimental and simulation data can be generated by drawing a sample from this jointly normal distribution.^{m} - (2)
Update the emulators after obtaining the hypothetical data.

To calculate how the emulators will change after adding the hypothetical data (which we treat as real data in the preposterior analysis), one could do another round of model bias correction using the collection of real data and hypothetical data, although this could be time consuming. Alternatively, by assuming that the estimation of the hyperparameters of the GP models will not significantly change with a few additional data, one can simply obtain a new model prediction by reusing the already estimated values of the hyperparameters in Eq. (5). The symbols now have new meanings; for example,

**d**and**V**in Eq. (5) are now referred to as the collection of real and hypothetical data and their covariance matrix, respectively._{d} - (3)
Conduct uncertainty propagation using the new model prediction.

After all the emulators for disciplinary responses are updated using the hypothetical data, uncertainty propagation as presented in Sec. 3.1 can provide information about how much the uncertainty of system QOIs

**y**_{sys}is changed. - (4)
Calculate the “expected” reduced uncertainty of system QOIs

**y**_{sys}.Note that steps 1–3 in this subsection are done under a specific set of hypothetical data. In order to achieve statistically meaningful results, steps 1–3 should be repeated by generating multiple sets of hypothetical data, which leads to an expected (averaged) reduced uncertainty of

**y**_{sys}.

The preposterior analysis can be conducted for any suggested resource allocation plan, and users can select a plan that balances effectiveness and cost. Practically, one can first analyze the expected uncertainty reduction under the plan that all additional resources are from simulations, and then check whether the criterion (9) will be satisfied. Those locations where this criterion is not likely to be satisfied after gathering new simulations should be considered for experimental measurements.

## Case Study

To demonstrate the proposed resource allocation approach, we consider an electronic packaging design problem [46], a benchmark multidisciplinary problem that has been frequently studied in the literature [23,25,47–49]. A circuit consisting of two resistors is mounted on a heat sink, and there are two coupled disciplines, i.e., the electrical and thermal subsystems, as demonstrated in Fig. 5. Eight input variables *x*_{1}–*x*_{8}, five linking variables *y*_{6}, *y*_{7}, and *y*_{11}–*y*_{13}, and a system QOI *y*_{1} are involved in this test problem. The numbering of the variables follows the same as used in the aforementioned works. Table 1 provides the physical meanings of the variables; a detailed problem statement that includes descriptions of all variables and the bounds of the design variables can be found in the referenced literature.

x_{1} | Heat sink width (m) |

x_{2} | Heat sink length (m) |

x_{3} | Fin length (m) |

x_{4} | Fin width (m) |

x_{5} | Nominal resistance #1 at temperature 20 °C (Ω) |

x_{6} | Temperature coefficient of electrical resistance #1 (K^{−1}) |

x_{7} | Nominal resistance #2 at temperature 20 °C (Ω) |

x_{8} | Temperature coefficient of electrical resistance #2 (K^{−1}) |

y_{1} | Negative of watt density (W/m^{3}) |

y_{6} | Power dissipation in resistor #1 (W) |

y_{7} | Power dissipation in resistor #2 (W) |

y_{11} | Component temperature of resistor #1 (°C) |

y_{12} | Component temperature of resistor #2 (°C) |

y_{13} | Heat sink volume (m^{3}) |

x_{1} | Heat sink width (m) |

x_{2} | Heat sink length (m) |

x_{3} | Fin length (m) |

x_{4} | Fin width (m) |

x_{5} | Nominal resistance #1 at temperature 20 °C (Ω) |

x_{6} | Temperature coefficient of electrical resistance #1 (K^{−1}) |

x_{7} | Nominal resistance #2 at temperature 20 °C (Ω) |

x_{8} | Temperature coefficient of electrical resistance #2 (K^{−1}) |

y_{1} | Negative of watt density (W/m^{3}) |

y_{6} | Power dissipation in resistor #1 (W) |

y_{7} | Power dissipation in resistor #2 (W) |

y_{11} | Component temperature of resistor #1 (°C) |

y_{12} | Component temperature of resistor #2 (°C) |

y_{13} | Heat sink volume (m^{3}) |

The original problem is a design optimization formulation to maximize the watt density (i.e., to minimize *y*_{1}). To test our proposed approach, we modify the original problem as follows. Instead of finding an optimal solution, we now aim at achieving a good predictive model of *y*_{1} over the whole design space of *x*_{1}–*x*_{8}. Since the thermal discipline is much more complex than the electrical discipline in that it requires a finite difference strategy to calculate the component temperatures, we assume that the models to evaluate *y*_{11}, *y*_{12}, given {*x*_{1}, *x*_{2}, *x*_{3}, *x*_{4}, *y*_{6}, *y*_{7}}, are inaccurate and require model bias correction. In the first iteration of resource allocation, we collect simulation and experimental data and conduct model bias correction for *y*_{11} and *y*_{12} (Sec. 4.1), and then evaluate the impact of epistemic model uncertainty on the system QOI *y*_{1} (Sec. 4.2), and then apply the proposed resource allocation approach to reduce the aggregated uncertainty on *y*_{1} (Sec. 4.3). The results from subsequent iterations are provided in Sec. 4.4.

### Model Bias Correction for the Thermal Discipline (First Iteration).

We first collect data comprised of 40 “experimental” observations. Since the original benchmark problem does not provide any experimental data, we generate a set of data for *y*_{11} and *y*_{12} by adding random noise to the exact model provided in Ref. [46], which we treat as the experimental data. The 40 experimental runs are determined via an LHS over the six-dimensional input space {*x*_{1}, *x*_{2}, *x*_{3}, *x*_{4}, *y*_{6}, *y*_{7}}. We then collect data comprised of 40 simulation observations of *y*_{11} and *y*_{12} from low-fidelity simulation models that we build, which intentionally have a significant discrepancy **δ**(**x**) compared with the exact model in Ref. [46]. The 40 simulation runs are determined via a different LHS, so that their locations are completely different than those for the experimental runs. Note that {*y*_{6}, *y*_{7}} are input variables to the thermal discipline and yet essentially linking variables in the whole system; hence, we do not know beforehand their ranges. We assume in this stage that they are both within [0, 8] based on some preliminary system-level simulations (using the current “wrong” simulation models). We will show later that their actual bounds are slightly beyond that.

The data points are shown in Fig. 6, with only the {*y*_{6}, *y*_{7}} plotted for illustration purpose. The GP-based model bias correction in Sec. 2.1 is applied to integrate the two sets of data and to subsequently obtain the updated emulators of *y*_{11} and *y*_{12.} An extra validation step is performed with 25 reserved experimental data; a comparison between the updated emulators and the validation data, together with the 95% PIs, is provided in Fig. 7. It can be seen that for both responses, 40 experimental data + 40 simulations are barely sufficient for bias correction: at some validation points, the PIs are large, and a few validation data are even outside the 95% PI. It can be anticipated that such uncertainty will propagate and influence the system QOI *y*_{1}.

### Selection of Input Settings and Responses for Resources Allocation (First Iteration).

In this step, 2000 LHS samples are generated to explore the input space of *x*_{1}–*x*_{8} (see Fig. 8). The MUA method described in Sec. 3.1 evaluates the aggregated uncertainty of *y*_{1} at each sample point. As a reminder, the analysis in this step only relies upon the emulators that we just built in Sec. 4.1, and thus, the computations for 2000 samples are not at all prohibitive.

Uncertainty acceptance criterion *α*% in Eq. (9) is set to be 10%. It turns out that there are only 11 samples (the highlighted blue rectangle and red circles in Fig. 8) out of 2000 that violate the acceptance criterion, which indicates an overall acceptable amount of uncertainty in the design space. However, the largest *γ* in Eq. (9) is 47.31%, which shows that in some local areas, the uncertainty is huge. Using the correlation check presented in Sec. 3.3 (with the correlation limit *β* set to be 0.95), 10 samples (the red circles) out of 11 are selected for allocating more resources.

The result of MSSA at the ten selected input settings is provided in Fig. 9, based on which we decide to allocate more resources to *y*_{11} at input settings #1–3, 5, and 10 (since at these points the MSI of *y*_{11} is larger than that of *y*_{12}), and to *y*_{12} at the rest of input settings. Note that here we only plot the MSI of *y*_{11} and *y*_{12}. It is because they are only two uncertainty sources in this problem, and based on the definitions of MSI and TSI in Eqs. (18) and (19), we have

We omit the plots of TSI as MSI offers sufficient information from sensitivity analysis.

### Preposterior Analysis to Decide the Type of Resources to Allocate (First Iteration).

We now apply the preposterior analysis to decide which type of resources to allocate, for each of the ten selected input settings. A total of 500 random sets of data are drawn from the joint normal distribution of $y11m$ at input settings #1–3, 5, and 10, and $y12m$ at input settings #4 and 6–9, which compose the hypothetical simulation data (a 500 × 10 matrix). This joint distribution comes from the GP models of $y11m$ and $y12m$ built in Sec. 4.1. By adding one set of hypothetical data (one row of the matrix), the emulators for *y*_{11} and *y*_{12} are both updated (hypothetically), based on which the aggregated uncertainty of system QOI *y*_{1} (presented by *γ* in Eq. (9)) is reevaluated (again, for each of the ten selected input settings). The predicted reduced uncertainty of *y*_{1} after gathering more simulation data is taken to be the average of the reevaluated *γ* over the 500 hypothetical sets of simulation data.

The result of the preposterior analysis is shown in Fig. 10. The uncertainty of *y*_{1} at most of the input settings (#2, 3, and 6–10) is expected to be reduced to the acceptable level (*α*% = 10% in Eq. (9)); however, points #1, 4, and 5 still do not satisfy the uncertainty criterion. Therefore, instead of allocating simulation data to all the ten points, we should rather allocate experimental data to #1, 4, and 5.

To sum up, our decision for resource allocation after the first iteration of analysis is:

- (1)
conduct simulations of

*y*_{11}at points #2, 3, and 10 and*y*_{12}at points #6–9 - (2)
conduct physical experiments of

*y*_{11}at points #1 and 5 and*y*_{12}at point #4

### Decisions for Resource Allocation (Second–Fourth Iterations).

The process of resource allocation is iterative. The 2000 LHS samples in each iteration are different, and hence help to thoroughly explore the whole design space. It turns out that after four iterations, all of the LHS samples generated satisfy the uncertainty limit specified in Eq. (9). The locations and types of the added resources and the corresponding responses after each iteration are provided in Fig. 11. The number of total samples (for model bias correction) after each iteration is provided in Table 2.

Iteration # | $y11e$ | $y11m$ | $y12e$ | $y12m$ |
---|---|---|---|---|

0 | 40 | 40 | 40 | 40 |

1 | 42 | 43 | 41 | 44 |

2 | 43 | 49 | 44 | 49 |

3 | 44 | 50 | 45 | 52 |

4 | 44 | 51 | 46 | 53 |

Iteration # | $y11e$ | $y11m$ | $y12e$ | $y12m$ |
---|---|---|---|---|

0 | 40 | 40 | 40 | 40 |

1 | 42 | 43 | 41 | 44 |

2 | 43 | 49 | 44 | 49 |

3 | 44 | 50 | 45 | 52 |

4 | 44 | 51 | 46 | 53 |

In Fig. 11, it is seen that although {*y*_{6}, *y*_{7}} are initially assumed to be within the range of [0, 8]^{2} (discussed in Sec. 4.1), the actual range for {*y*_{6}, *y*_{7}} is slightly beyond that and there are some new samples selected in the later iterations that reach the region where *y*_{6} > 8 or *y*_{7} > 8. We consider it a merit of the proposed method that it can gradually detect the region where we might have overlooked as the models get improved.

Beyond considering only physical experiments in existing literature, we treat computer simulations as yet another useful tool to reduce the epistemic model uncertainty, and this example showcases the merit of allocating computational resources. As shown in Table 2, more than 70% of the allocated resources (24 out of 34) are from computer simulations. We believe that it is because computer simulations help reduce the uncertainty due to the lack of data, and potentially also help discover the fundamental features of disciplinary responses. One can imagine that in an ideal situation when simulation models are perfect, uncertainty reduction can be achieved by conducting only simulations without help from physical experiments.

Another observation from this case study is that, to achieve a globally satisfactory model for QOIs in a multidisciplinary system does not require the disciplinary subsystem emulators to be globally satisfactory in their own input domains. In other words, the allocated resources to disciplinary subsystems do not necessarily cover the input space of disciplinary models. For example, it can be seen from Fig. 11 that all the resources allocated to *y*_{11} are clustered in a narrow band where *y*_{7} is small, and that the resources allocated to *y*_{12} are clustered in a narrow band where *y*_{6} is small. It is due to the complex coupling relationship between different disciplines in the system; many combinations of {*y*_{6}, *y*_{7}} are not achievable when the two disciplines couple together, and hence, there is no need to improve the modeling capability over those unachievable regions.

## Conclusions

In this paper, we develop a new approach for resource allocation in multidisciplinary systems to reduce the aggregated epistemic model uncertainty of system QOIs. The novelty of the approach lies in breaking down the complex resource allocation problem that involves multiple decisions into a sequential set of decisions to answer questions of *where* (sampling locations), *what* (disciplinary responses), and *which* (simulations versus experiments). The proposed sequential decision-making strategy manages the complexity in creating globally accurate surrogate models for each subsystem in multidisciplinary analysis and design optimization by allocating resources (samples) that have the largest impact on reducing the uncertainty in predicting system QOI. After quantifying disciplinary epistemic model uncertainty using GP-based model bias correction, we first explore the input space of the whole system to identify the locations (*where*) with unacceptable amounts of uncertainty with respect to the system QOIs. In this step, the MUA method provides a fast and analytical uncertainty analysis. Next, the identification of critical responses (*what*) at the selected inputs is done using an MSSA we developed earlier to evaluate the relative contribution of each functional response. Meanwhile, the proposed correlation analysis uses the information from GP-based correlation to determine the proximity of input settings and to ensure that the selected input settings are sparsely located. Finally, decisions are made about *which* type of resources to allocate to the critical responses at the chosen input locations, via an enhanced preposterior analysis that predicts the effect of gathering more resources prior to the actual resource allocation. The preposterior analysis, unlike any preceding work in the literature, is applicable to any combination of experiments and simulations. The proposed method is applied to a benchmark electronic packaging problem to reduce the aggregated epistemic model uncertainty of the system QOI.

The proposed approach has the following advantages. First, it considers not only the physical experiments but also the computer simulations. We demonstrate in the case study that in many cases, the physical experiments are not really necessary, if the epistemic model uncertainty mainly comes from lack of data. Second, because the proposed approach strategically breaks resource allocation into a sequential process, the decision making is much more tractable. Third, the method is efficient for complex multidisciplinary analysis. Except for the stage when the simulation or physical experiments are collected, all the analyses are conducted based on inexpensive GP emulators. The MUA and MSSA methods further accelerate the procedure by providing analytical formulas. Last but not least, while we do not elaborate it, the proposed framework can be simplified to work for a single discipline as well. In that case, MUA and MSSA may not be useful as a single discipline typically involves a single response; however, the GP-based correlation analysis and the preposterior analysis are still applicable.

The limitation of the work is that it does not explicitly consider the costs of simulations and experiments, although users can adjust the tunable preference parameters of the uncertainty acceptance level and correlation coefficient of samples based on the total budget. Our future research will include the budget as a constraint in the decision-making process of resource allocation.

## Acknowledgment

The grant support from the National Science Foundation (CMMI-1233403 and CMMI-1537641) is greatly acknowledged. Shishi Chen's predoctoral visit at Northwestern University is sponsored by the China Scholarship Council.

## Nomenclature

**d**=collected data

**h**(^{m}**·**),**h**(^{δ}**·**) =predefined basis function for polynomial regression of mean function of simulation model GP/discrepancy function GP

- MSI =
main sensitivity index

- ND =
number of disciplines

- TSI =
total sensitivity index

**u**=_{ij}linking variables, output from the

*i*th discipline and input to the*j*th discipline*V*(·,·),^{m}*V*(·,·) =^{δ}covariance function of simulation model GP/discrepancy function GP

**x**=(

*x*_{1},…,*x*) =_{p}*p*-dimensional input variable**x**,_{i}**x**_{ind}=disciplinary input variables of the

*i*th discipline/all ND disciplines**x**=_{s}shared input variables across two or more disciplines

**y**,_{i}**y**_{ind}=disciplinary outputs of the

*i*th discipline/all ND disciplines**y**_{sys}=system QOIs

*y*(·),^{m}*y*(·) =^{e}simulation/experimental response

- $y\u0302e$ (
**·**) =mean prediction of the experimental response

*Z*(**·**) =a random quantity representing model uncertainty

*Z*=_{l}model uncertainty term for model

*L***Z**=_{∼l}model uncertainty terms for models except

*L**α*% =predefined uncertainty limit for

*γ*(·,·)*β*=predefined limit for

*ρ***β**,^{m}**β**=^{δ}coefficients for polynomial regression of mean function of simulation model GP/discrepancy function GP

*γ*(·,·) =measure of uncertainty at a particular input setting

*δ*(·) =model discrepancy function

*ε*=experimental error

**μ**=mean vector; subscript, if available, denotes a specific random quantify, e.g.,

**μ**. refers to the mean of_{ui}**u**._{i}*ρ*=correlation coefficient

- $\sigma Z2$ (
**·**) =variance of

*Z*(**·**) **Σ**=covariance matrix; subscript, if available, denotes a specific random quantify, e.g.,

**Σ**refers to the covariance_{u}**u**^{e}**ϕ**=hyperparameters of GP

### Appendix

##### Pseudocode for Choosing the Sample Points in Sec. 3.3.

Given: $(xind(k),xs(k))$, k = 1,…, N from Section 3.1 _{s} | |

1. | Sort k by $Var[ysys(xind(k),xs(k))]$ in descending order |

2. | Initialize α, β, list1 = ${\xi |(xind(\xi ),xs(\xi ))violatesEq.(9)}$, list2 = {1} |

3. | for all ξ∈ list1 |

4. | for all η ∈ list2 |

5. | find response L that contributes most to $Var[ysys(xind(\xi ),xs(\xi ))]$ |

6. | find response L′ that contributes most to $Var[ysys(xind(\eta ),xs(\eta ))]$ |

7. | if L ≠ L′, then |

8. | set flag(ξ,η) = 1 |

9. | elseif ρ_{ξ}_{,} < _{η}β [Eq. (20)], then |

10. | set flag(ξ,η) = 1 |

11. | else |

12. | set flag(ξ,η) = 0 |

13. | endif |

14. | endfor |

15. | if flag(ξ,η) = 1, ∀η∈ list2, then |

16. | add ξ to list2 |

17. | endif |

18. | endfor |

Given: $(xind(k),xs(k))$, k = 1,…, N from Section 3.1 _{s} | |

1. | Sort k by $Var[ysys(xind(k),xs(k))]$ in descending order |

2. | Initialize α, β, list1 = ${\xi |(xind(\xi ),xs(\xi ))violatesEq.(9)}$, list2 = {1} |

3. | for all ξ∈ list1 |

4. | for all η ∈ list2 |

5. | find response L that contributes most to $Var[ysys(xind(\xi ),xs(\xi ))]$ |

6. | find response L′ that contributes most to $Var[ysys(xind(\eta ),xs(\eta ))]$ |

7. | if L ≠ L′, then |

8. | set flag(ξ,η) = 1 |

9. | elseif ρ_{ξ}_{,} < _{η}β [Eq. (20)], then |

10. | set flag(ξ,η) = 1 |

11. | else |

12. | set flag(ξ,η) = 0 |

13. | endif |

14. | endfor |

15. | if flag(ξ,η) = 1, ∀η∈ list2, then |

16. | add ξ to list2 |

17. | endif |

18. | endfor |

The chosen sample points are stored in “list2.”