Model distance-based approach for global sensitivity analysis in engineering systems with multivariate outputs

Abstract

In engineering systems, uncertainties in input parameters can significantly influence the output responses. This paper proposes a model distance-based approach to perform global sensitivity analysis for quantifying the influence of input uncertainties on multiple responses in an engineering system. The sensitivity indices are determined by comparing a reference model that incorporates all system uncertainties, with an altered model, where specific uncertainties are constrained. The proposed framework employs probability distance measures such as Hellinger distance, Kullback–Leibler divergence, and l₂ norm which are based on joint probability density functions. The study also demonstrates the equivalence between the l₂ norm-based approach and Sobol’s analysis in multivariate sensitivity context. The proposed methodology effectively manages correlated random variables, accommodates both Gaussian and non-Gaussian distributions, and allows for the grouping of input variables. Illustrative examples consist of static analysis of a truss system and dynamic analysis of a frame subjected to seismic excitation. The sensitivity indices are estimated using brute-force Monte Carlo simulations. The relative ranking of these sensitivity indices can be utilized to identify the most and least significant variables contributing to the response uncertainty. The numerical results show a consistent ranking of input variables across different probability measures, indicating the robustness of proposed framework.

Graphical abstract

Keywords

multivariate outputs / moment-independent sensitivity analysis / probability distance measures / multivariate Sobol’s analysis / Monte Carlo simulation / uncertainty quantification

1 Introduction

Uncertainty is an inherent feature of engineering systems, which arises from incomplete knowledge of input parameters, unpredictable environmental conditions, and complex interactions within the system. Fundamentally, uncertainty results from the absence of precise information about the behavior of the system, which can be better understood through fundamental concepts such as experiments, trials, and outcomes. An experiment refers to any physical phenomenon that can be repeatedly observed, with each repetition constituting a trial, and each trial producing an outcome. While the range of possible outcomes for an experiment may be known, the specific outcome of any individual trial remains uncertain. This uncertainty in trial outcomes leads to uncertainty in the observed variables. Uncertainties in engineering systems include parameter, model, and data uncertainties, which are inherent to many engineering problems. As a result, forecasting how this system will respond, its cost, and its performance under different circumstances is more difficult. However, these diverse uncertainties can be categorized as either aleatory or epistemic uncertainty [1]. Probabilistic theory can be used to model aleatory uncertainty, whereas non-probabilistic theories such as interval, convex, or fuzzy concepts can be used to model epistemic uncertainty [2]. The response of the mathematical or computational model is influenced by the uncertainty involved in the input parameters. The level of uncertainty in the response of the model is quantified using sensitivity analysis. This helps in using the engineering systems efficiently while admitting its limitations and associated uncertainty [3]. In this regard, Saltelli et al. [4] sees that sensitivity analysis approaches used in the context of uncertainty analysis need to conform to the underlying three criteria: “quantitative, global, and model-free.” The term global refers to the ability of the method to account for the overall input distribution. Model independence shows that the sensitivity analysis technique can give reliable outcomes without making any assumptions about how the model is dependent on its input parameters. Many fields, such as reliability analysis, reliability optimization, and reliability design, use sensitivity analysis [5].

Local sensitivity analysis and global sensitivity analysis are two primary categories of methods used in sensitivity analysis. With a given input space, local sensitivity analysis helps to determine which input parameters have a significant effect on the response and how sensitive the response is to small changes in those parameters [4]. The partial derivatives of the response with respect to the inputs serve as the foundation for these techniques. The assessment of uncertainty in the output is limited to a certain reference value. This entails systematically adjusting input random variables individually, one at a time, and by small magnitudes. This analysis loses its potential to offer global sensitivity information [6]. In contrast, global sensitivity analysis methods quantify the impact of uncertainties in the input random variable across the entire input space of the model. This analysis is an effective tool for identifying significant input parameters, evaluating the relationships among input parameters, assessing the robustness of the model and uncertainty, and improving comprehension of the model [7]. Global sensitivity analysis has the potential to improve the dependability, accuracy, and decision-making capacity of engineering systems by achieving these goals [3].

The various methods to perform global sensitivity analysis include the screening method [8], elementary effect method [9], variance-based method [10–13], and moment-independent method [14]. The variance-based method and moment-independent method are widely used approaches in global sensitivity analysis. Sobol’s method is a variance-based method that provides a decomposition of the response variance into contributions from individual input parameters and their interactions with other input parameters [15]. Specifically, Sobol indices quantify the proportion of variance attributed to single or multiple input variables. These indices serve as quantitative measures of the influence exerted by input variables on the variance of response of the model [16]. This technique is known as analysis of variance. It was initially developed for scalar response with input parameters that are uniformly distributed and follow a sequence of independent and identically distributed random variables [7,15]. This notion has been expanded and generalized over time to deal with different conditions, such as dependent inputs [17–19] and so on. Sobol’s method has been extensively studied and has found application in modeling of engineering systems [20], often regarded as an exact solution and commonly used to validate the accuracy of new methodologies [21]. Within the context of structural engineering, this method is extensively examined across static and dynamic as well as linear/nonlinear problem frameworks [15,22–24]. Kucherenko et al. [25] suggested using derivative-based method to reduce computational costs resulting from the higher dimensionality of the input vector in global sensitivity analysis. This sensitivity measure is performed by averaging the partial derivatives over the input variable space.

Global sensitivity measures in moment-independent methods account for the complete distribution of the model response, such as the cumulative distribution function (CDF) or the probability density function (PDF), rather than focusing on specific moments. This approach has been widely employed in various engineering scenarios including chemical reactor studies [26], investigations of biological systems [27], environmental modeling [28], and hydrological modeling [29]. When dealing with skewed or multi-modal probability density function of the response, these strategies are especially useful as they effectively capture the decision-maker’s level of knowledge about the model response [20]. These moment-independent techniques can be performed using some of the probability distance measures. Chun et al. [30] gave the Minkowski distance metric based on the CDFs of the model response. Greegar and Manohar [31] employed probability measures such as Hellinger distance, Kullback–Leibler divergence, and l₂ norm to perform sensitivity analysis on scalar responses. The authors also showed that l₂ norm-based indices are related to Sobol’s indices. Greegar and Manohar [32] also broadened the scope of global sensitivity analysis based on model distance, in the context of non-probabilistic methodologies for modeling uncertain input variables. Nandi and Singh [33] used probability measures, such as the Kantorovich–Rubinstein metric, Hellinger distance, total variation distance, Kolmogorov, Bhattacharya, and Cramer-von Mises for scalar response. Borgonovo [14] identified a new measure called the moment-independent uncertainty indicator (delta index) which studies the changes in the density functions attributed due to the changes in the uncertainty characteristics of the input variables. Abhinav and Manohar [34] extended moment-independent approach in modeling the vibrating systems subjected to random excitations. They introduced model distance-based sensitivity measures concerning the response power spectral density function and PDF. Greegar et al. [35] extended the studies on model distance-based global-local sensitivity analysis to linear structural dynamic problems with time-dependent random excitations and spatially varying stochastic inhomogeneities. Li et al. [36] proposed PDF-based sensitivity measure to analyze the uncertainties in seismic demand of bridges.

In the context of models with multivariate outputs, traditional global sensitivity analysis on individual responses or specific scalar functions often falls short [37]. Multivariate global sensitivity analysis addresses this limitation by evaluating how input parameters simultaneously influence multiple and interdependent responses. This analysis captures interactions between responses that single-variate methods ignore, this prevents incomplete or misleading conclusions about which input parameters are most influential. It reduces the computational complexity of conducting multiple independent analyses for each output by offering a unified framework to analyze input parameter effects across all responses. Also, multivariate global sensitivity analysis facilitates the identification and ranking of critical input variables by providing a more precise and comprehensive evaluation of their influence on multiple system outputs. This approach enhances the accuracy of sensitivity assessments, which enables a deeper understanding of the contributions of each input to the overall system behavior. Zhang et al. [38] used the variance decomposition method for performing a global sensitivity analysis of multivariate outputs. There are many methods available for global sensitivity analysis of multivariate outputs, such as multivariate analysis of variance [39], covariance-based sensitivity analysis [40], principal component analysis [41,42], distance components decomposition [43], and distance correlation-based approach [44]. Gamboa et al. [40] proposed the covariance decomposition method for multivariate outputs, which includes decomposing the covariance of model response. The proposed sensitivity index in this method was later shown to be identical to those defined in the variance-based decomposition method. Li et al. [45] presented a sensitivity index based on probability integral transformation, which involves transforming any continuous distribution of a random variable to a standard uniform distribution using CDF. This method assesses the significance of each input variable by examining its influence on the joint CDFs of response of the model. Energy distance-based sensitivity analysis also utilizes joint CDFs [46,47]. Liu et al. [48] introduced a multiple response Gaussian process surrogate model with separable covariance to estimate sensitivity indices for multiple responses. Sun et al. [49] suggested a method to efficiently estimate projection-based multivariate sensitivity indices using two polynomial chaos-based surrogates: polynomial chaos expansion and a proper orthogonal decomposition-based polynomial chaos expansion. Liu et al. [50] proposed a double-loop Kriging method to estimate multivariate sensitivity indices. Chen and Huang [51] proposed a generalized hybrid metamodel combining radial basis function and polynomial chaos expansion to accommodate various distribution types for multivariate outputs. Tang et al. [52] developed a multi-output global sensitivity analysis method based on a correlation-analysis-guided Kriging model that addresses high-dimensional inputs and multivariate outputs in high-speed train dynamics. Chiani et al. [53] introduced a global sensitivity analysis method for integrated assessment models with multivariate outputs using optimal transport theory. Mazo and Tournier [54] developed a multivariate global sensitivity analysis method based on statistical inference on total indices using Möbius transform techniques.

The present study was motivated by the following observations.

1) Model distance-based global sensitivity analysis for scalar response using probability distance measures provides an opportunity to extend this concept to multivariate response.

2) For scalar output, when the input random variables are independent, l₂ norm-based indices are related to Sobol’s indices. Since the Sobol multivariate global sensitivity analysis relies on the second-order moment [40], it prompts an idea to investigate the equivalence between sensitivity indices obtained from Sobol’s analysis and l₂ norm-based approach.

3) Model distance-based multivariate sensitivity analysis can be used for evaluating sensitivity measures with respect to grouped random variables.

The first objective of this paper is to perform multivariate global sensitivity analysis based on model distance-based approach by utilizing joint PDFs of the responses. These joint PDFs are estimated using multivariate kernel density estimator and brute-force Monte Carlo simulation (MCS). Probability distance measures such as Hellinger distance, Kullback–Leibler divergence, and l₂ norm are used for this evaluation. This analysis is conducted across various scenarios, including uncorrelated, correlated, and grouped input random variables. The second objective of this paper is to show the equivalence between l₂ norm-based method and Sobol’s method in the context of multivariate global sensitivity indices, when the input random variables are independent. To demonstrate these proposed concepts, the study considers a truss carrying uncertain loads with uncertain parameters and a two-story steel frame subjected to ground motion. The rest of this paper is structured as follows. Section 2 provides an overview of Sobol’s analysis for multivariate outputs. Section 3 discusses the model distance-based approach for multivariate sensitivity analysis. Section 4 shows the equivalence between l₂ norm-based approach and Sobol’s method. Section 5 discusses the model distance-based sensitivity analysis for grouped variables. A brief procedure for evaluation of global sensitivity indices using brute-force MCS is discussed in Section 6. Section 7 contains the demonstration of the proposed sensitivity indices, and Section 8 highlights the significant findings of the proposed framework, followed by the conclusions.

2 Multivariate Sobol’s analysis

This section provides an overview of Sobol’s analysis for multivariate response. Let \mathit{X} represent n-dimensional input random vector (i.e., \mathit{X}={\left[X_1,X_2,\dots ,X_n\right]}^{\mathrm{T}}) with marginal PDFs, p_{X_i}\left(x_i\right);i=1,2,\dots ,n, and joint PDF, p_{\mathit{X}}\left(\mathit{x}\right)=\prod _{i=1}^n{p}_{X_i}\left(x_i\right). The superscript T represents the transpose of a vector. Consider a deterministic function Y^{\left(r\right)}=f^{\left(r\right)}\left(\mathit{X}\right);r=1,2,\dots ,s where the function f^{\left(r\right)}(\cdot ) could be the result of a simulation of a mathematical model and Y^{\left(\mathrm{r}\right)} denotes the response of the model. Let \mathit{Y} be s-dimensional response vector (i.e., \mathit{Y}=[Y^{\left(1\right)},Y^{\left(2\right)},\dots ,Y^{\left(s\right)}{]}^{\mathrm{t}}) with joint PDF, p_{\mathit{Y}}\left(\mathit{y}\right). Sobol’s decomposition of Y^{\left(r\right)} is presented as follows [40].

(1)

with

(2)f_0^{\left(r\right)}=\mathrm{E}\left[Y^{\left(r\right)}\right]=\underset{-\mathrm{\infty }}{\stackrel{\mathrm{\infty }}{\int }}\cdots \underset{-\mathrm{\infty }}{\stackrel{\mathrm{\infty }}{\int }}{f}^{\left(r\right)}\left(X_1,X_2,\dots ,X_n\right)\prod _{k=1}^n{p}_{X_k}\left(x_k\right)\mathrm{d}{x}_k,

(3)\begin{array}{rl}{f}_i^{\left(r\right)}\left(X_i\right)=& \mathrm{E}\left[Y^{\left(r\right)}|X_i\right]-\mathrm{E}\left[Y^{\left(r\right)}\right]\\ =& \underset{-\mathrm{\infty }}{\stackrel{\mathrm{\infty }}{\int }}\cdots \underset{-\mathrm{\infty }}{\stackrel{\mathrm{\infty }}{\int }}{f}^{\left(r\right)}\left(X_1,X_2,\dots ,X_n\right)\prod _{\begin{array}{c}k=1\\ k\ne i\end{array}}^n{p}_{X_k}\left(x_k\right)\mathrm{d}{x}_k-f_0^{\left(r\right)},\end{array}

(4)\begin{array}{rl}{f}_{ij}^{\left(r\right)}\left(X_i,X_j\right)=& \mathrm{E}\left[Y^{\left(r\right)}|X_i,X_j\right]-\mathrm{E}\left[Y^{\left(r\right)}|X_i\right]-\mathrm{E}\left[Y^{\left(r\right)}|X_j\right]-\mathrm{E}\left[Y^{\left(r\right)}\right]\\ =& \underset{-\mathrm{\infty }}{\stackrel{\mathrm{\infty }}{\int }}\cdots \underset{-\mathrm{\infty }}{\stackrel{\mathrm{\infty }}{\int }}{f}^{\left(r\right)}\left(X_1,X_2,\dots ,X_n\right)\prod _{\begin{array}{c}k=1\\ k\ne i\\ k\ne j\end{array}}^n{p}_{X_k}\left(x_k\right)\mathrm{d}{x}_k\\ & -f_0^{\left(r\right)}-f_i^{\left(r\right)}\left(X_i\right)-f_j^{\left(r\right)}\left(X_j\right),\end{array}

(5)\begin{array}{rl}{f}_{12\cdots n}^{\left(r\right)}\left(X_1,X_2,\dots ,X_n\right)=& f^{\left(r\right)}\left(X_1,X_2,\dots ,X_n\right)\\ & -\left(\mathrm{s}\mathrm{u}\mathrm{m}\mathrm{m}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{s}\mathrm{o}\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{e}\mathrm{r}\mathrm{o}\mathrm{r}\mathrm{d}\mathrm{e}\mathrm{r}\right),\end{array}

where \mathrm{E}[\cdot ] denotes the mathematical expectation operator. The terms in Eq. (1) should satisfy the orthogonal properties as follows [31].

(6)\underset{-\mathrm{\infty }}{\stackrel{\mathrm{\infty }}{\int }}\cdots \underset{-\mathrm{\infty }}{\stackrel{\mathrm{\infty }}{\int }}{f}_{i_1{i}_2\cdots i_m}^{\left(r\right)}\left(x_{i_1},x_{i_2},\dots ,x_{i_m}\right)\prod _{k=1}^n{p}_{X_k}\left(x_k\right)\mathrm{d}{x}_k=0,

(7)\begin{array}{rl}& \underset{-\mathrm{\infty }}{\stackrel{\mathrm{\infty }}{\int }}\cdots \underset{-\mathrm{\infty }}{\stackrel{\mathrm{\infty }}{\int }}{f}_{i_1{i}_2\cdots i_m}^{\left(r\right)}\left(x_{i_1},x_{i_2},\dots ,x_{i_m}\right)f_{j_1{j}_2\cdots j_q}^{\left(r\right)}\left(x_{j_1},x_{j_2},\dots ,x_{j_q}\right)\\ & \prod _{k=1}^n{p}_{X_k}\left(x_k\right)\mathrm{d}{x}_{\mathrm{k}}\{\begin{array}{l}\ne 0,\mathrm{i}\mathrm{f}m=q\mathrm{a}\mathrm{n}\mathrm{d}{i}_{\mathrm{l}}=j_{\mathrm{l}}\mathrm{\forall }l=1,2,\dots ,n,\\ =0,\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e},\end{array}\end{array}

(8)\underset{-\mathrm{\infty }}{\stackrel{\mathrm{\infty }}{\int }}{f}_{i_1{i}_2\cdots i_m}^{\left(r\right)}\left(x_{i_1},x_{i_2},\dots ,x_{i_m}\right)p_{X_k}\left(x_k\right)\mathrm{d}{x}_k=0;\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{y}1\le k\le n.

For Eqs. (6)−(8), the condition is that ; 1\le m\le n. In Eq. (8), the integral of f_{i_1{i}_2\cdots i_m}^{\left(r\right)}(x_{i_1},x_{i_2},\dots ,x_{i_m}) with respect to any of its variables is equal to zero. Using Eq. (1), Gamboa et al. [40] represented Sobol’s decomposition of \mathit{Y} as

(9)

where

(10)\begin{array}{rl}& {\mathit{f}}_0=\left\{\begin{array}{c}\begin{array}{c}{f}_0^{\left(1\right)}\\ f_0^{\left(2\right)}\end{array}\\ \cdots \\ f_0^{\left(s\right)}\end{array}\right\},{\mathit{Y}}_{X_i}=\left[\begin{array}{c}{f}_i^{\left(1\right)}\left(X_i\right)\\ f_i^{\left(2\right)}\left(X_i\right)\\ \begin{array}{c}\cdots \\ f_i^{\left(s\right)}\left(X_i\right)\end{array}\end{array}\right],{\mathit{Y}}_{X_i{X}_j}=\left[\begin{array}{c}{f}_{ij}^{\left(1\right)}\left(X_i,X_j\right)\\ f_{ij}^{\left(2\right)}\left(X_i,X_j\right)\\ \begin{array}{c}\cdots \\ f_{ij}^{\left(s\right)}\left(X_i,X_j\right)\end{array}\end{array}\right],\\ & \mathrm{a}\mathrm{n}\mathrm{d}{\mathit{Y}}_{X_1{X}_2\cdots X_n}=\left[\begin{array}{c}{f}_{12\cdots n}^{\left(1\right)}\left(X_1,X_2,\dots ,X_n\right)\\ f_{12\cdots n}^{\left(2\right)}\left(X_1,X_2,\dots ,X_n\right)\\ \begin{array}{c}⋮\\ f_{12\cdots n}^{\left(s\right)}\left(X_1,X_2,\dots ,X_n\right)\end{array}\end{array}\right].\end{array}

By taking covariance on both sides of Eq. (9) gives

(11)

where

(12)\begin{array}{rl}\mathrm{C}\mathrm{o}\mathrm{v}\left(\mathit{Y}\right)=& \mathrm{E}\left[\left(\mathit{Y}-\mathrm{E}\left[\mathit{Y}\right]\right){\left(\mathit{Y}-\mathrm{E}\left[\mathit{Y}\right]\right)}^{\mathrm{t}}\right]\\ =& \mathrm{E}\left[\left\{\begin{array}{c}\begin{array}{c}{Y}^{\left(1\right)}-\mathrm{E}\left[Y^{\left(1\right)}\right]\\ Y^{\left(2\right)}-\mathrm{E}\left[Y^{\left(2\right)}\right]\end{array}\\ \cdots \\ Y^{\left(s\right)}-\mathrm{E}\left[Y^{\left(s\right)}\right]\end{array}\right\}{\left\{\begin{array}{c}\begin{array}{c}{Y}^{\left(1\right)}-\mathrm{E}\left[Y^{\left(1\right)}\right]\\ Y^{\left(2\right)}-\mathrm{E}\left[Y^{\left(2\right)}\right]\end{array}\\ \cdots \\ Y^{\left(s\right)}-\mathrm{E}\left[Y^{\left(s\right)}\right]\end{array}\right\}}^{\mathrm{t}}\right],\end{array}

where \mathrm{C}\mathrm{o}\mathrm{v}\left(\mathit{Y}\right) represents the covariance matrix of \mathit{Y}. In general, {\mathrm{C}\mathrm{o}\mathrm{v}}_{X_i}(\cdot ), {\mathrm{C}\mathrm{o}\mathrm{v}}_{X_i{X}_j}(\cdot ), and {\mathrm{C}\mathrm{o}\mathrm{v}}_{X_1{X}_2\cdots X_n}(\cdot ) can be written as follows

(13)\begin{array}{rl}{\mathrm{C}\mathrm{o}\mathrm{v}}_{X_i}\left({\mathit{Y}}_{X_i}\right)=& \mathrm{E}\left[\left({\mathit{Y}}_{X_i}-\mathrm{E}\left[{\mathit{Y}}_{X_i}\right]\right){\left({\mathit{Y}}_{X_i}-\mathrm{E}\left[{\mathit{Y}}_{X_i}\right]\right)}^{\mathrm{t}}\right]\\ =& \left[\begin{array}{llll}\mathrm{C}\mathrm{o}\mathrm{v}\left(f_i^{\left(1\right)}\left(X_i\right),f_i^{\left(1\right)}\left(X_i\right)\right)& \mathrm{C}\mathrm{o}\mathrm{v}\left(f_i^{\left(1\right)}\left(X_i\right),f_i^{\left(2\right)}\left(X_i\right)\right)& \cdots & \mathrm{C}\mathrm{o}\mathrm{v}\left(f_i^{\left(1\right)}\left(X_i\right),f_i^{\left(s\right)}\left(X_i\right)\right)\\ \mathrm{C}\mathrm{o}\mathrm{v}\left(f_i^{\left(2\right)}\left(X_i\right),f_i^{\left(1\right)}\left(X_i\right)\right)& \mathrm{C}\mathrm{o}\mathrm{v}\left(f_i^{\left(2\right)}\left(X_i\right),f_i^{\left(2\right)}\left(X_i\right)\right)& \cdots & \mathrm{C}\mathrm{o}\mathrm{v}\left(f_i^{\left(2\right)}\left(X_i\right),f_i^{\left(s\right)}\left(X_i\right)\right)\\ \cdots & & & \\ \mathrm{C}\mathrm{o}\mathrm{v}\left(f_i^{\left(s\right)}\left(X_i\right),f_i^{\left(1\right)}\left(X_i\right)\right)& \mathrm{C}\mathrm{o}\mathrm{v}\left(f_i^{\left(s\right)}\left(X_i\right),f_i^{\left(2\right)}\left(X_i\right)\right)& \cdots & \mathrm{C}\mathrm{o}\mathrm{v}\left(f_i^{\left(s\right)}\left(X_i\right),f_i^{\left(s\right)}\left(X_i\right)\right)\end{array}\right],\end{array}

(14){\mathrm{C}\mathrm{o}\mathrm{v}}_{X_i{X}_j}\left({\mathit{Y}}_{X_i{X}_j}\right)=\left[\begin{array}{llll}\mathrm{C}\mathrm{o}\mathrm{v}\left(f_{ij}^{\left(1\right)}\left(X_i,X_j\right),f_{ij}^{\left(1\right)}\left(X_i,X_j\right)\right)& \mathrm{C}\mathrm{o}\mathrm{v}\left(f_{ij}^{\left(1\right)}\left(X_i,X_j\right),f_{ij}^{\left(2\right)}\left(X_i,X_j\right)\right)& \cdots & \mathrm{C}\mathrm{o}\mathrm{v}\left(f_{ij}^{\left(1\right)}\left(X_i,X_j\right),f_{ij}^{\left(s\right)}\left(X_i,X_j\right)\right)\\ \mathrm{C}\mathrm{o}\mathrm{v}\left(f_{ij}^{\left(2\right)}\left(X_i,X_j\right),f_{ij}^{\left(1\right)}\left(X_i,X_j\right)\right)& \mathrm{C}\mathrm{o}\mathrm{v}\left(f_{ij}^{\left(2\right)}\left(X_i,X_j\right),f_{ij}^{\left(2\right)}\left(X_i,X_j\right)\right)& \cdots & \mathrm{C}\mathrm{o}\mathrm{v}\left(f_{ij}^{\left(2\right)}\left(X_i,X_j\right),f_{ij}^{\left(s\right)}\left(X_i,X_j\right)\right)\\ \cdots & & & \\ \mathrm{C}\mathrm{o}\mathrm{v}\left(f_{ij}^{\left(s\right)}\left(X_i,X_j\right),f_{ij}^{\left(1\right)}\left(X_i,X_j\right)\right)& \mathrm{C}\mathrm{o}\mathrm{v}\left(f_{ij}^{\left(s\right)}\left(X_i,X_j\right),f_{ij}^{\left(2\right)}\left(X_i,X_j\right)\right)& \cdots & \mathrm{C}\mathrm{o}\mathrm{v}\left(f_{ij}^{\left(s\right)}\left(X_i,X_j\right),f_{ij}^{\left(s\right)}\left(X_i,X_j\right)\right)\end{array}\right],

(15)\begin{array}{rl}& {\mathrm{C}\mathrm{o}\mathrm{v}}_{X_1{X}_2\cdots X_n}\left({\mathit{Y}}_{X_1{X}_2\cdots X_n}\right)=\\ & \left[\begin{array}{llll}\mathrm{C}\mathrm{o}\mathrm{v}\left(f_{12\cdots n}^{\left(1\right)}\left(X_1,X_2,\dots ,X_n\right),f_{12\cdots n}^{\left(1\right)}\left(X_1,X_2,\dots ,X_n\right)\right)& \mathrm{C}\mathrm{o}\mathrm{v}\left(f_{12\cdots n}^{\left(1\right)}\left(X_1,X_2,\dots ,X_n\right),f_{12\cdots n}^{\left(2\right)}\left(X_1,X_2,\dots ,X_n\right)\right)& \cdots & \mathrm{C}\mathrm{o}\mathrm{v}\left(f_{12\cdots n}^{\left(1\right)}\left(X_1,X_2,\dots ,X_n\right),f_{12\cdots n}^{\left(s\right)}\left(X_1,X_2,\dots ,X_n\right)\right)\\ \mathrm{C}\mathrm{o}\mathrm{v}\left(f_{12\cdots n}^{\left(2\right)}\left(X_1,X_2,\dots ,X_n\right),f_{12\cdots n}^{\left(1\right)}\left(X_1,X_2,\dots ,X_n\right)\right)& \mathrm{C}\mathrm{o}\mathrm{v}\left(f_{12\cdots n}^{\left(2\right)}\left(X_1,X_2,\dots ,X_n\right),f_{12\cdots n}^{\left(2\right)}\left(X_1,X_2,\dots ,X_n\right)\right)& \cdots & \mathrm{C}\mathrm{o}\mathrm{v}\left(f_{12\cdots n}^{\left(2\right)}\left(X_1,X_2,\dots ,X_n\right),f_{12\cdots n}^{\left(s\right)}\left(X_1,X_2,\dots ,X_n\right)\right)\\ \cdots & & & \\ \mathrm{C}\mathrm{o}\mathrm{v}\left(f_{12\cdots n}^{\left(s\right)}\left(X_1,X_2,\dots ,X_n\right),f_{12\cdots n}^{\left(1\right)}\left(X_1,X_2,\dots ,X_n\right)\right)& \mathrm{C}\mathrm{o}\mathrm{v}\left(f_{12\cdots n}^{\left(s\right)}\left(X_1,X_2,\dots ,X_n\right),f_{12\cdots n}^{\left(2\right)}\left(X_1,X_2,\dots ,X_n\right)\right)& \cdots & \mathrm{C}\mathrm{o}\mathrm{v}\left(f_{12\cdots n}^{\left(s\right)}\left(X_1,X_2,\dots ,X_n\right),f_{12\cdots n}^{\left(s\right)}\left(X_1,X_2,\dots ,X_n\right)\right)\end{array}\right],\end{array}

where the term {\mathrm{C}\mathrm{o}\mathrm{v}}_{X_i}(\cdot ) represents the covariance matrix of \mathit{Y} solely influenced by input variable X_i. The terms {\mathrm{C}\mathrm{o}\mathrm{v}}_{X_i{X}_j}(\cdot ) and {\mathrm{C}\mathrm{o}\mathrm{v}}_{X_1{X}_2\cdots X_n}(\cdot ) denote the covariance matrix of \mathit{Y} influenced by the pair (X_i,X_j) and nth order tuples, respectively. By taking the trace on both sides of Eq. (11) gives

(16)

Equation (16) can be rewritten as

(17)

Equation (17) can be denoted as

(18a)

(18b)\begin{array}{c}{D}^{\mathit{Y}}=\mathrm{V}\mathrm{a}\mathrm{r}\left[Y^{\left(1\right)}\right]+\mathrm{V}\mathrm{a}\mathrm{r}\left[Y^{\left(2\right)}\right]+\cdots +\mathrm{V}\mathrm{a}\mathrm{r}\left[Y^{\left(s\right)}\right],\end{array}

(18c)\begin{array}{rl}\sum _{i=1}^n{D}_i=& \sum _{i=1}^n\{\mathrm{V}\mathrm{a}\mathrm{r}\left[f_i^{\left(1\right)}\left(X_i\right)\right]+\mathrm{V}\mathrm{a}\mathrm{r}\left[f_i^{\left(2\right)}\left(X_i\right)\right]+\cdots \\ & +\mathrm{V}\mathrm{a}\mathrm{r}\left[f_i^{\left(s\right)}\left(X_i\right)\right]\},\end{array}

(18d)

(18e)\begin{array}{rl}{D}_{12\cdots n}=& \{\mathrm{V}\mathrm{a}\mathrm{r}\left[f_{12\cdots n}^{\left(1\right)}\left(X_1,X_2,\dots ,X_n\right)\right]\\ & +\mathrm{V}\mathrm{a}\mathrm{r}\left[f_{12\cdots n}^{\left(2\right)}\left(X_1,X_2,\dots ,X_n\right)\right]+\cdots \\ & +\mathrm{V}\mathrm{a}\mathrm{r}\left[f_{12\cdots n}^{\left(s\right)}\left(X_1,X_2,\dots ,X_n\right)\right]\},\end{array}

where D^{\mathit{Y}} denotes the variance of \mathit{Y}. The term D_{\mathrm{i}} represents the partial variance of \mathit{Y} contributed by input X_i alone. The terms D_{ij} and D_{12\cdots n} denote the partial variance of \mathit{Y} contributed by the pair (X_i,X_j) and nth order tuples, respectively. Using Eq. (18a), Sobol’s first order and total sensitivity indices of input variable X_i, which are denoted as S_i and T_i respectively, can be defined as follows [40].

(19)S_i=\frac{D_i}{D^{\mathit{Y}}};i=1,2,\dots ,n,

(20)T_i=\frac{D_i+\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^n{D}_{ij}+\cdots +D_{12\cdots n}}{D^{\mathit{Y}}};i=1,2,\dots ,n,

where S_i indicates how each input variable impacts the overall variability of \mathit{Y} by assessing its contribution to the overall response variance and T_i represents the variable’s overall impact on overall variability of \mathit{Y}, by considering both the variable and its interaction with other input variables. Higher values of S_i and T_i imply higher sensitivity effect of variable X_i toward the uncertainty in the response quantity of interest.

2.1 Remarks

1) If T_i=0, then it indicates that the ith input variable, either directly or by its interactions with other variables, is not contributing to the variability of the response \mathit{Y}.

2) If T_i=S_i=1, then it implies that \mathit{Y} is solely dependent on X_i only.

3) If S_i=0, then it implies that X_i is not directly contributing to the variability of the response \mathit{Y}.

3 Model distance-based sensitivity analysis for multivariate response

This section proposes model distance-based sensitivity analysis for multivariate response and shows the equivalence between probability distance measure based on l₂ norm and Sobol’s method in the context of multivariate global sensitivity analysis. Consider the reference model, \mathit{Y}=f\left(\mathit{X}\right). All the input variables of \mathit{X} in this model are treated as uncertain. Consider another hypothetical model, in which all the input parameters are treated as uncertain except the ith input parameter, X_i which has been deliberately fixed at a deterministic value denoted by {\alpha }_i. Let this model be called altered model-1. The response of this altered model-1 is denoted by {\mathit{Y}}_i. Let p_{\mathit{Y}}\left(y^{\left(1\right)},y^{\left(2\right)},\dots ,y^{\left(s\right)}\right) denotes the joint PDF of Y and p_{{\mathit{Y}}_i}\left(y_i^{\left(1\right)},y_i^{\left(2\right)},\dots ,y_i^{\left(s\right)}\right) denotes the joint PDF of {\mathit{Y}}_i. The difference between the reference model, \mathit{Y} and the altered model-1, {\mathit{Y}}_i represents the importance of the eliminated uncertainty in X_i and its interactions with other input variables. Let this difference be denoted by d_i=\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{t}(\mathit{Y},{\mathit{Y}}_i) and can be considered as a measure of total sensitivity index of the input variable, X_i.

Similarly, an altered model-2 can be defined where except for the ith element, X_i, all other variables are made deterministic by fixing them to constant values, {\alpha }_j;j=1,2,\dots ,n;j\ne i. Fig.1 illustrates the reference model, and the altered models 1 and 2. Let the response of this altered model-2 is denoted by {\tilde{\mathit{Y}}}_i and the corresponding joint PDF is given by p_{{\tilde{\mathit{Y}}}_i}\left({\tilde{y}}_i^{\left(1\right)},{\tilde{y}}_i^{\left(2\right)},\dots ,{\tilde{y}}_i^{\left(s\right)}\right). The difference between the reference model, Y and the altered model-2, {\tilde{\mathit{Y}}}_i represents the sensitivity of the uncertainty due to X_i alone. This sensitivity measure denoted by d_{\sim i}=\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{t}\left(\mathit{Y},{\tilde{\mathit{Y}}}_i\right) can be considered as the main sensitivity index due to variable X_i.

The difference between the reference model and the hypothetical models can be measured by comparing the joint PDFs/CDFs using various probability distance measures. In the existing literature, different probability distance measures are available to evaluate the disparity between two random variables by comparing their PDFs/CDFs [20,55], such as Hellinger distance, Kullback–Leibler divergence, and l₂ norm. This study focuses on extending these probability measures aimed at measuring the distance between two random vectors. Tab.1 illustrates a few probability distance measures used to compare random vectors.

In Tab.1, U and V represent the random vectors with joint PDFs, denoted by p_{\mathit{U}}(u_1,u_2,\dots ,u_s) and p_{\mathit{V}}(u_1,u_2,\dots ,u_s), respectively, and with joint CDFs, denoted by P_{\mathit{U}}(u_1,u_2,\dots ,u_s) and P_{\mathit{V}}(u_1,u_2,\dots ,u_s), respectively.

The probability measures listed in Tab.1 can be used for defining the total and main sensitivity indices by comparing the joint PDFs of Y, {\mathit{Y}}_i, and {\tilde{\mathit{Y}}}_i. For instance, to find the distance between Y and {\mathit{Y}}_i, using Hellinger distance, Kullback–Leibler divergence, and l₂ norm are defined in Eqs. (21)−(23).

(21)d_{\mathrm{H}}\left(\mathit{Y},{\mathit{Y}}_i\right)=\sqrt{\frac{1}{2}\underset{-\mathrm{\infty }}{\stackrel{\mathrm{\infty }}{\int }}\cdots \underset{-\mathrm{\infty }}{\stackrel{\mathrm{\infty }}{\int }}{\left(\sqrt{p_{\mathbf{Y}}\left(y^{\left(1\right)},y^{\left(2\right)},\dots ,y^{\left(s\right)}\right)}-\sqrt{p_{{\mathbf{Y}}_i}\left(y^{\left(1\right)},y^{\left(2\right)},\dots ,y^{\left(s\right)}\right)}\right)}^2\prod _{i=1}^s\mathrm{d}{y}^{\left(i\right)}},

(22)d_{\mathrm{K}}\left(\mathit{Y},{\mathit{Y}}_i\right)=\underset{-\mathrm{\infty }}{\stackrel{\mathrm{\infty }}{\int }}\cdots \underset{-\mathrm{\infty }}{\stackrel{\mathrm{\infty }}{\int }}{p}_{\mathbf{Y}}\left(y^{\left(1\right)},y^{\left(2\right)},\dots ,y^{\left(s\right)}\right)\mathrm{l}\mathrm{n}\frac{p_{\mathbf{Y}}\left(y^{\left(1\right)},y^{\left(2\right)},\dots ,y^{\left(s\right)}\right)}{p_{{\mathbf{Y}}_i}\left(y^{\left(1\right)},y^{\left(2\right)},\dots ,y^{\left(s\right)}\right)}\prod _{i=1}^s\mathrm{d}{y}^{\left(i\right)},

(23)l_2\left(\mathit{Y},{\mathit{Y}}_i\right)=\sqrt{\mathrm{t}\mathrm{r}\left[\mathrm{E}\left[\left(\mathit{Y}-{\mathit{Y}}_i\right){\left(\mathit{Y}-{\mathit{Y}}_i\right)}^{\mathrm{t}}\right]\right]},

where Eqs. (21)−(23) represent the total sensitivity indices evaluated with respect to the input random variable, X_i. Similarly, the main sensitivity indices with respect to the input random variable, X_i based on Hellinger distance, Kullback–Leibler divergence, and l₂ norm, can also be defined by using d_{\mathrm{H}}\left(\mathit{Y},{\tilde{\mathit{Y}}}_{\mathrm{i}}\right), d_{\mathrm{K}}\left(\mathit{Y},{\tilde{\mathit{Y}}}_{\mathrm{i}}\right), and l₂\left(\mathit{Y},{\tilde{\mathit{Y}}}_{\mathrm{i}}\right), respectively.

3.1 Remarks

1) Any distance measure is said to be a metric if it satisfies the properties of non-negativity, identity, symmetry, and triangle inequality. Here, Hellinger distance and l₂ norm satisfy all these properties, and can be considered as a metric. However, Kullback–Leibler divergence, while satisfying the non-negativity and identity properties, does not satisfy symmetry property and triangle inequality. Hence, it is considered as a divergence measure not as a metric.

2) The d_i calculated with respect to input variable X_i represents the degree of uncertainty that X_i imparts along with its interactions with other input variables to the response. A higher value of d_i suggests that the variable X_i has more significant sensitivity effect.

3) The d_{\sim i} evaluated with respect to input variable X_i represents the level of uncertainty X_i alone holds toward the response. Higher the value of d_{\sim i} indicates that uncertainty in X_i has less importance on the response.

4) The probability distance measures, d_i and d_{\sim i} evaluated using Hellinger distance, fall within the range of 0 to 1. The distance measures evaluated using Kullback–Leibler divergence are greater than or equal to 0.

4 Equivalence between l₂ norm and Sobol’s multivariate global sensitivity analysis

Sobol’s multivariate global sensitivity analysis relies on the second-order moment (variance) of the response variables [40]. In such a case, it becomes relevant to investigate the equivalence between the resulting sensitivity indices from Sobol’s analysis and the metric, l₂ norm (i.e. {\left\{l_2\left(\mathit{Y},{\mathit{Y}}_i\right)\right\}}^2=\mathrm{E}\left[\left(\mathit{Y}-{\mathit{Y}}_i\right){\left(\mathit{Y}-{\mathit{Y}}_i\right)}^{\mathrm{t}}\right]) in which the sensitivity indices also depend only on the second-order moment of the response variables. To illustrate this, consider a deterministic function, Y^{\left(r\right)}=f^{\left(r\right)}\left(X_1,X_2,X_3\right); r=1\mathrm{a}\mathrm{n}\mathrm{d}2, where \mathit{Y}={\left[Y^{\left(1\right)},Y^{\left(2\right)}\right]}^{\mathrm{t}} with joint PDF p_{\mathit{Y}}\left(y_1,y_2\right) and X_1, X_2, and X_3 are independent random variables. The joint PDF of \mathit{X} is denoted as p_{\mathit{X}}\left(x_1,x_2,x_3\right) which is the product of marginal PDFs of \mathit{X}, \prod _{i=1}^3{p}_{X_i}\left(x_i\right). Sobol’s decomposition of Y^{\left(1\right)} is represented as follows

(24)\begin{array}{rl}{Y}^{\left(1\right)}=& f_0^{\left(1\right)}+f_1^{\left(1\right)}\left(X_1\right)+f_2^{\left(1\right)}\left(X_2\right)+f_3^{\left(1\right)}\left(X_3\right)+f_{12}^{\left(1\right)}\left(X_1,X_2\right)\\ & +f_{13}^{\left(1\right)}\left(X_1,X_3\right)+f_{23}^{\left(1\right)}\left(X_2,X_3\right)+f_{123}^{\left(1\right)}\left(X_1,X_2,X_3\right).\end{array}

The summands of Sobol’s decomposition of Y^{\left(1\right)} are evaluated as follows

(25)f_0^{\left(1\right)}=\underset{-\mathrm{\infty }}{\stackrel{\mathrm{\infty }}{\int }}\underset{-\mathrm{\infty }}{\stackrel{\mathrm{\infty }}{\int }}\underset{-\mathrm{\infty }}{\stackrel{\mathrm{\infty }}{\int }}{f}^{\left(1\right)}\left(X_1,X_2,X_3\right)\prod _{k=1}^3{p}_{X_k}\left(x_k\right)\mathrm{d}{x}_k,

(26)\begin{array}{rl}{f}_i^{\left(1\right)}\left(X_i\right)=& \underset{-\mathrm{\infty }}{\stackrel{\mathrm{\infty }}{\int }}\underset{-\mathrm{\infty }}{\stackrel{\mathrm{\infty }}{\int }}\underset{-\mathrm{\infty }}{\stackrel{\mathrm{\infty }}{\int }}{f}^{\left(1\right)}\left(X_1,X_2,X_3\right)\\ & \cdot \prod _{\begin{array}{c}k=1\\ k\ne i\end{array}}^3{p}_{X_k}\left(x_k\right)\mathrm{d}{x}_k-f_0^{\left(1\right)};i=1,2,3,\end{array}

(27)

(28)\begin{array}{rl}{f}_{123}^{\left(1\right)}\left(X_1,X_2,X_3\right)=& f^{\left(1\right)}\left(X_1,X_2,X_3\right)\\ & -\left(\mathrm{s}\mathrm{u}\mathrm{m}\mathrm{m}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{s}\mathrm{o}\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{e}\mathrm{r}\mathrm{o}\mathrm{r}\mathrm{d}\mathrm{e}\mathrm{r}\right).\end{array}

In a similar manner, the summands of Sobol’s decomposition of Y^{\left(2\right)}can also be evaluated. Then, using Eq. (9), Sobol’s decomposition of \mathit{Y} can be represented as follows

(29)

Following the notation used in Eq. (10) with n = 3, the above equation can be rewritten as

(30)\begin{array}{rl}\left\{\begin{array}{c}{Y}^{\left(1\right)}\\ Y^{\left(2\right)}\end{array}\right\}=& \left\{\begin{array}{c}{f}_0^{\left(1\right)}\\ f_0^{\left(2\right)}\end{array}\right\}+\left\{\begin{array}{c}{f}_1^{\left(1\right)}\left(X_1\right)\\ f_1^{\left(2\right)}\left(X_1\right)\end{array}\right\}+\left\{\begin{array}{c}{f}_2^{\left(1\right)}\left(X_2\right)\\ f_2^{\left(2\right)}\left(X_2\right)\end{array}\right\}\\ & +\left\{\begin{array}{c}{f}_3^{\left(1\right)}\left(X_3\right)\\ f_3^{\left(2\right)}\left(X_3\right)\end{array}\right\}+\left\{\begin{array}{c}{f}_{12}^{\left(1\right)}\left(X_1,X_2\right)\\ f_{12}^{\left(2\right)}\left(X_1,X_2\right)\end{array}\right\}+\left\{\begin{array}{c}{f}_{13}^{\left(1\right)}\left(X_1,X_3\right)\\ f_{13}^{\left(2\right)}\left(X_1,X_3\right)\end{array}\right\}\\ & +\left\{\begin{array}{c}{f}_{23}^{\left(1\right)}\left(X_2,X_3\right)\\ f_{23}^{\left(2\right)}\left(X_2,X_3\right)\end{array}\right\}+\left\{\begin{array}{c}{f}_{123}^{\left(1\right)}\left(X_1,X_2,X_3\right)\\ f_{123}^{\left(2\right)}\left(X_1,X_2,X_3\right)\end{array}\right\}.\end{array}

Taking covariance on both sides of Eq. (29) gives

(31)

where

(32)\begin{array}{rl}\mathrm{C}\mathrm{o}\mathrm{v}\left(\mathit{Y}\right)=& \mathrm{E}\left[\left\{\begin{array}{c}{Y}^{\left(1\right)}-\mathrm{E}\left[Y^{\left(1\right)}\right]\\ Y^{\left(2\right)}-\mathrm{E}\left[Y^{\left(2\right)}\right]\end{array}\right\}{\left\{\begin{array}{c}{Y}^{\left(1\right)}-\mathrm{E}\left[Y^{\left(1\right)}\right]\\ Y^{\left(2\right)}-\mathrm{E}\left[Y^{\left(2\right)}\right]\end{array}\right\}}^{\mathrm{t}}\right]\\ =& \left[\begin{array}{cc}\mathrm{V}\mathrm{a}\mathrm{r}\left[Y^{\left(1\right)}\right]& \mathrm{C}\mathrm{o}\mathrm{v}\left(Y^{\left(1\right)},Y^{\left(2\right)}\right)\\ \mathrm{C}\mathrm{o}\mathrm{v}\left(Y^{\left(1\right)},Y^{\left(2\right)}\right)& \mathrm{V}\mathrm{a}\mathrm{r}\left[Y^{\left(2\right)}\right]\end{array}\right].\end{array}

In general, {\mathrm{C}\mathrm{o}\mathrm{v}}_{X_i}\left({\mathit{Y}}_{X_i}\right) and {\mathrm{C}\mathrm{o}\mathrm{v}}_{X_i{X}_j}\left({\mathit{Y}}_{X_i{X}_j}\right) can be written as given by Eqs. (33) and (34).

(33){\mathrm{C}\mathrm{o}\mathrm{v}}_{X_i}\left({\mathit{Y}}_{X_i}\right)=\left[\begin{array}{cc}\mathrm{V}\mathrm{a}\mathrm{r}\left[f_i^{\left(1\right)}\left(X_i\right)\right]& \mathrm{C}\mathrm{o}\mathrm{v}\left(f_i^{\left(1\right)}\left(X_i\right),f_i^{\left(2\right)}\left({\mathrm{X}}_{\mathrm{i}}\right)\right)\\ \mathrm{C}\mathrm{o}\mathrm{v}\left(f_i^{\left(1\right)}\left(X_i\right),f_i^{\left(2\right)}\left(X_i\right)\right)& \mathrm{V}\mathrm{a}\mathrm{r}\left[f_i^{\left(2\right)}\left({\mathrm{X}}_{\mathrm{i}}\right)\right]\end{array}\right];i=1,2,\mathrm{a}\mathrm{n}\mathrm{d}3,

(34)

(35)\mathrm{C}\mathrm{o}{\mathrm{v}}_{X_1{X}_2{X}_3}\left({\mathit{Y}}_{X_1{X}_2{X}_3}\right)=\left[\begin{array}{cc}\mathrm{V}\mathrm{a}\mathrm{r}\left[f_{123}^{\left(1\right)}\left(X_1,X_2,X_3\right)\right]& \mathrm{C}\mathrm{o}\mathrm{v}\left(f_{123}^{\left(1\right)}\left(X_1,X_2,X_3\right),f_{123}^{\left(2\right)}\left(X_1,X_2,X_3\right)\right)\\ \mathrm{C}\mathrm{o}\mathrm{v}\left(f_{123}^{\left(1\right)}\left(X_1,X_2,X_3\right),f_{123}^{\left(2\right)}\left(X_1,X_2,X_3\right)\right)& \mathrm{V}\mathrm{a}\mathrm{r}\left[f_{123}^{\left(2\right)}\left(X_1,X_2,X_3\right)\right]\end{array}\right].

Taking trace on both sides of Eq. (31) gives

(36)

This further reduces as follows

(37)\begin{array}{rl}\mathrm{V}\mathrm{a}\mathrm{r}\left[Y^{\left(1\right)}\right]& +\mathrm{V}\mathrm{a}\mathrm{r}\left[Y^{\left(2\right)}\right]=\mathrm{V}\mathrm{a}\mathrm{r}\left[f_1^{\left(1\right)}\left(X_1\right)\right]+\mathrm{V}\mathrm{a}\mathrm{r}\left[f_1^{\left(2\right)}\left(X_1\right)\right]\\ & +\mathrm{V}\mathrm{a}\mathrm{r}\left[f_2^{\left(1\right)}\left(X_2\right)\right]+\mathrm{V}\mathrm{a}\mathrm{r}\left[f_2^{\left(2\right)}\left(X_2\right)\right]\\ & +\mathrm{V}\mathrm{a}\mathrm{r}\left[f_3^{\left(1\right)}\left(X_3\right)\right]+\mathrm{V}\mathrm{a}\mathrm{r}\left[f_3^{\left(2\right)}\left(X_3\right)\right]\\ & +\mathrm{V}\mathrm{a}\mathrm{r}\left[f_{12}^{\left(1\right)}\left(X_1,X_2\right)\right]+\mathrm{V}\mathrm{a}\mathrm{r}\left[f_{12}^{\left(2\right)}\left(X_1,X_2\right)\right]\\ & +\mathrm{V}\mathrm{a}\mathrm{r}\left[f_{13}^{\left(1\right)}\left(X_1,X_3\right)\right]+\mathrm{V}\mathrm{a}\mathrm{r}\left[f_{13}^{\left(2\right)}\left(X_1,X_3\right)\right]\\ & +\mathrm{V}\mathrm{a}\mathrm{r}\left[f_{23}^{\left(1\right)}\left(X_2,X_3\right)\right]+\mathrm{V}\mathrm{a}\mathrm{r}\left[f_{23}^{\left(2\right)}\left(X_2,X_3\right)\right]\\ & +\mathrm{V}\mathrm{a}\mathrm{r}\left[f_{123}^{\left(1\right)}\left(X_1,X_2,X_3\right)\right]+\mathrm{V}\mathrm{a}\mathrm{r}\left[f_{123}^{\left(2\right)}\left(X_1,X_2,X_3\right)\right].\end{array}

Therefore, the total variance in Y can be written as

(38)D^{\mathit{Y}}=D_1+D_2+D_3+D_{12}+D_{13}+D_{23}+D_{123},

where D_1=\mathrm{V}\mathrm{a}\mathrm{r}\left[f_1^{\left(1\right)}\left(X_1\right)\right]+\mathrm{V}\mathrm{a}\mathrm{r}\left[f_1^{\left(2\right)}\left(X_1\right)\right], which represents the sum of partial variances of responses Y^{\left(1\right)}\mathrm{a}\mathrm{n}\mathrm{d}{Y}^{\left(2\right)} with respect to the input variable X_1 alone. In a similar manner, D_2\mathrm{a}\mathrm{n}\mathrm{d}{D}_3 can be evaluated which denote the sum of partial variances of responses Y^{\left(1\right)}\mathrm{a}\mathrm{n}\mathrm{d}{Y}^{\left(2\right)} with respect to the variables, X_2{\mathrm{a}\mathrm{n}\mathrm{d}X}_3 respectively. The higher order partial variances of responses Y^{\left(1\right)}\mathrm{a}\mathrm{n}\mathrm{d}{Y}^{\left(2\right)} are denoted by D_{12},D_{13},D_{23}, and D_{123}. Using Eqs. (19) and (20), Sobol’s first order and total sensitivity indices of X_i; i=1,2,3 can be quantified as follows

(39)S_1=\frac{D_1}{D^{\mathit{Y}}},S_2=\frac{D_2}{D^{\mathit{Y}}},\mathrm{a}\mathrm{n}\mathrm{d}{S}_3=\frac{D_3}{D^{\mathit{Y}}}.

(40)\begin{array}{rl}& T_1=\frac{D_1+D_{12}+D_{13}+D_{123}}{D^{\mathit{Y}}},T_2=\frac{D_2+D_{12}+D_{23}+D_{123}}{D^{\mathit{Y}}},\\ & \mathrm{a}\mathrm{n}\mathrm{d}{T}_3=\frac{D_3+D_{13}+D_{23}+D_{123}}{D^{\mathit{Y}}}.\end{array}

To determine the total sensitivity index based on l₂ norm, begin by fixing X_1 at a constant value \alpha. Using Eq. (29), {\mathit{Y}}_1 can be represented as follows

(41)\begin{array}{rl}{\mathit{Y}}_1=& \left\{\begin{array}{c}{Y}_1^{\left(1\right)}\\ Y_1^{\left(2\right)}\end{array}\right\}=\left\{\begin{array}{c}{f}_0^{\left(1\right)}\\ f_0^{\left(2\right)}\end{array}\right\}+\left\{\begin{array}{c}{f}_1^{\left(1\right)}\left(\alpha \right)\\ f_1^{\left(2\right)}\left(\alpha \right)\end{array}\right\}+\left\{\begin{array}{c}{f}_2^{\left(1\right)}\left(X_2\right)\\ f_2^{\left(2\right)}\left(X_2\right)\end{array}\right\}+\\ & \left\{\begin{array}{c}{f}_3^{\left(1\right)}\left(X_3\right)\\ f_3^{\left(2\right)}\left(X_3\right)\end{array}\right\}+\left\{\begin{array}{c}{f}_{12}^{\left(1\right)}\left(\alpha ,X_2\right)\\ f_{12}^{\left(2\right)}\left(\alpha ,X_2\right)\end{array}\right\}+\left\{\begin{array}{c}{f}_{13}^{\left(1\right)}\left(\alpha ,X_3\right)\\ f_{13}^{\left(2\right)}\left(\alpha ,X_3\right)\end{array}\right\}\\ & +\left\{\begin{array}{c}{f}_{23}^{\left(1\right)}\left(X_2,X_3\right)\\ f_{23}^{\left(2\right)}\left(X_2,X_3\right)\end{array}\right\}+\left\{\begin{array}{c}{f}_{123}^{\left(1\right)}\left(\alpha ,X_2,X_3\right)\\ f_{123}^{\left(2\right)}\left(\alpha ,X_2,X_3\right)\end{array}\right\}.\end{array}

The l_2 norm between \mathit{Y} and {\mathit{Y}}_1 can be evaluated as follows

(42)\begin{array}{rl}{\left\{l_2\left(\mathit{Y},{\mathit{Y}}_1\right)\right\}}^2=& \mathrm{E}\left[\left(\mathit{Y}-{\mathit{Y}}_1\right){\left(\mathit{Y}-{\mathit{Y}}_1\right)}^{\mathrm{t}}\right]\\ =& \mathrm{E}\left[\left\{\begin{array}{c}{Y}^{\left(1\right)}-Y_1^{\left(1\right)}\\ Y^{\left(2\right)}-Y_1^{\left(2\right)}\end{array}\right\}{\left\{\begin{array}{c}{Y}^{\left(1\right)}-Y_1^{\left(1\right)}\\ Y^{\left(2\right)}-Y_1^{\left(2\right)}\end{array}\right\}}^{\mathrm{t}}\right],\end{array}

Taking trace of Eq. (42) gives

(43)\begin{array}{rl}& \mathrm{t}\mathrm{r}\left[{\left\{l_2\left(\mathit{Y},{\mathit{Y}}_1\right)\right\}}^2\right]=\mathrm{E}\left[{\left(Y^{\left(1\right)}-Y_1^{\left(1\right)}\right)}^2\right]+\mathrm{E}\left[{\left(Y^{\left(2\right)}-Y_1^{\left(2\right)}\right)}^2\right],\\ =& \mathrm{E}\left[{\left\{f_1^{\left(1\right)}\left(X_1\right)+f_{12}^{\left(1\right)}\left(X_1,X_2\right)+f_{13}^{\left(1\right)}\left(X_1,X_3\right)+f_{123}^{\left(1\right)}\left(X_1,X_2,X_3\right)-f_1^{\left(1\right)}\left(\alpha \right)-f_{12}^{\left(1\right)}\left(\alpha ,X_2\right)-f_{13}^{\left(1\right)}\left(\alpha ,X_3\right)-f_{123}^{\left(1\right)}\left(\alpha ,X_2,X_3\right)\right\}}^2\right]\\ & +\mathrm{E}\left[{\left\{f_1^{\left(2\right)}\left(X_1\right)+f_{12}^{\left(2\right)}\left(X_1,X_2\right)+f_{13}^{\left(2\right)}\left(X_1,X_3\right)+f_{123}^{\left(2\right)}\left(X_1,X_2,X_3\right)-f_1^{\left(2\right)}\left(\alpha \right)-f_{12}^{\left(2\right)}\left(\alpha ,X_2\right)-f_{13}^{\left(2\right)}\left(\alpha ,X_3\right)-f_{123}^{\left(2\right)}\left(\alpha ,X_2,X_3\right)\right\}}^2\right],\end{array}

where \mathrm{t}\mathrm{r}\left[{\left\{l_2\left(\mathit{Y},{\mathit{Y}}_1\right)\right\}}^2\right] can be rewritten as follows

(44)\begin{array}{rl}\mathrm{t}\mathrm{r}\left[{\left\{l_2\left(\mathit{Y},{\mathit{Y}}_1\right)\right\}}^2\right]=& D_1+D_{12}+D_{13}+D_{123}+{\left\{f_1^{\left(1\right)}\left(\alpha \right)\right\}}^2+{\left\{f_1^{\left(2\right)}\left(\alpha \right)\right\}}^2+\mathrm{E}\left[{\left\{f_{12}^{\left(1\right)}\left(\alpha ,X_2\right)\right\}}^2\right]+\mathrm{E}\left[{\left\{f_{12}^{\left(2\right)}\left(\alpha ,X_2\right)\right\}}^2\right]+\mathrm{E}\left[{\left\{f_{13}^{\left(1\right)}\left(\alpha ,X_3\right)\right\}}^2\right]\\ & +\mathrm{E}\left[{\left\{f_{13}^{\left(2\right)}\left(\alpha ,X_3\right)\right\}}^2\right]+\mathrm{E}\left[{\left\{f_{123}^{\left(1\right)}\left(\alpha ,X_2,X_3\right)\right\}}^2\right]+\mathrm{E}\left[{\left\{f_{123}^{\left(2\right)}\left(\alpha ,X_2,X_3\right)\right\}}^2\right].\end{array}

To formulate global total sensitivity index for input variable X_1, treat \alpha as a random variable, following the distributional properties of the random variable, X_1. Then evaluate \mathrm{E}\left[\mathrm{t}\mathrm{r}\left[{\left\{l_2\left(\mathit{Y},{\mathit{Y}}_1\right)\right\}}^2\right]\right] as follows

(45)\mathrm{E}\left[\mathrm{t}\mathrm{r}\left[{\left\{l_2\left(\mathit{Y},{\mathit{Y}}_1\right)\right\}}^2\right]\right]=2\left(D_1+D_{12}+D_{13}+D_{123}\right).

Using Eq. (40), the above expression can be written as

(46)\mathrm{E}\left[\mathrm{t}\mathrm{r}\left[{\left\{l_2\left(\mathit{Y},{\mathit{Y}}_1\right)\right\}}^2\right]\right]=2{D^{\mathit{Y}}T}_1.

Equation (46) represents the relationship between total sensitivity index with respect to the input variable X_1 based on l₂ norm and Sobol’s analysis. In a similar manner, total sensitivity indices with respect to input variables X_2 and X_3 can be evaluated based on l₂ norm. The general relationship between total sensitivity indices based on l₂ norm and Sobol’s method is shown in Eq. (47).

(47)T_i=\frac{\mathrm{E}\left[\mathrm{t}\mathrm{r}\left[{\left\{l_2\left(\mathit{Y},{\mathit{Y}}_i\right)\right\}}^2\right]\right]}{2D^{\mathit{Y}}};i=1,2,3.

Furthermore, by fixing X_2\mathrm{a}\mathrm{s}\alpha and X_3\mathrm{a}\mathrm{s}\beta, the main sensitivity index of X_1 can be evaluated using l₂ norm. Using Eq. (29), {\tilde{\mathit{Y}}}_1 can be represented as follows

(48)\begin{array}{rl}{\tilde{\mathit{Y}}}_1=& \left\{\begin{array}{c}{\tilde{Y}}_1^{\left(1\right)}\\ {\tilde{Y}}_1^{\left(2\right)}\end{array}\right\}=\left\{\begin{array}{c}{f}_0^{\left(1\right)}\\ f_0^{\left(2\right)}\end{array}\right\}+\left\{\begin{array}{c}{f}_1^{\left(1\right)}\left(X_1\right)\\ f_1^{\left(2\right)}\left(X_1\right)\end{array}\right\}+\left\{\begin{array}{c}{f}_2^{\left(1\right)}\left(\alpha \right)\\ f_2^{\left(2\right)}\left(\alpha \right)\end{array}\right\}\\ & +\left\{\begin{array}{c}{f}_3^{\left(1\right)}\left(\beta \right)\\ f_3^{\left(2\right)}\left(\beta \right)\end{array}\right\}+\left\{\begin{array}{c}{f}_{12}^{\left(1\right)}\left(X_1,\alpha \right)\\ f_{12}^{\left(2\right)}\left(X_1,\alpha \right)\end{array}\right\}+\left\{\begin{array}{c}{f}_{13}^{\left(1\right)}\left(X_1,\beta \right)\\ f_{13}^{\left(2\right)}\left(X_1,\beta \right)\end{array}\right\}\\ & +\left\{\begin{array}{c}{f}_{23}^{\left(1\right)}\left(\alpha ,\beta \right)\\ f_{23}^{\left(2\right)}\left(\alpha ,\beta \right)\end{array}\right\}+\left\{\begin{array}{c}{f}_{123}^{\left(1\right)}\left(X_1,\alpha ,\beta \right)\\ f_{123}^{\left(2\right)}\left(X_1,\alpha ,\beta \right)\end{array}\right\}.\end{array}

The l_2 norm between \mathit{Y} and {\tilde{\mathit{Y}}}_1 can be evaluated as follows

(49)\begin{array}{rl}{\left\{l_2\left(\mathit{Y},{\tilde{\mathit{Y}}}_1\right)\right\}}^2=& \mathrm{E}\left[\left(\mathit{Y}-{\tilde{\mathit{Y}}}_1\right){\left(\mathit{Y}-{\tilde{\mathit{Y}}}_1\right)}^{\mathrm{t}}\right],\\ =& \mathrm{E}\left[\left\{\begin{array}{c}{Y}^{\left(1\right)}-{\tilde{Y}}_1^{\left(1\right)}\\ Y^{\left(2\right)}-{\tilde{Y}}_1^{\left(2\right)}\end{array}\right\}\left\{\left(Y^{\left(1\right)}-{\tilde{Y}}_1^{\left(1\right)}\right)\left(Y^{\left(2\right)}-{\tilde{Y}}_1^{\left(2\right)}\right)\right\}\right].\end{array}

Taking trace of Eq. (49) gives

(50)\begin{array}{rl}\mathrm{t}\mathrm{r}\left[{\left\{l_2\left(\mathit{Y},{\tilde{\mathit{Y}}}_1\right)\right\}}^2\right]=& \mathrm{E}\left[{\left(Y^{\left(1\right)}-{\tilde{Y}}_1^{\left(1\right)}\right)}^2\right]+\mathrm{E}\left[{\left(Y^{\left(2\right)}-{\tilde{Y}}_1^{\left(2\right)}\right)}^2\right]=\mathrm{E}[\{f_2^{\left(1\right)}\left(X_2\right)+f_3^{\left(1\right)}\left(X_3\right)+f_{12}^{\left(1\right)}\left(X_1,X_2\right)+f_{13}^{\left(1\right)}\left(X_1,X_3\right)+f_{23}^{\left(1\right)}\left(X_2,X_3\right)\\ & {+f_{123}^{\left(1\right)}\left(X_1,X_2,X_3\right)-f_2^{\left(1\right)}\left(\alpha \right)-f_3^{\left(1\right)}\left(\beta \right)-f_{12}^{\left(1\right)}\left(X_1,\alpha \right)-f_{13}^{\left(1\right)}\left(X_1,\beta \right)-f_{23}^{\left(1\right)}\left(\alpha ,\beta \right)-f_{123}^{\left(1\right)}\left(X_1,\alpha ,\beta \right)\}}^2]\\ & +\mathrm{E}[\{f_2^{\left(2\right)}\left(X_2\right)+f_3^{\left(2\right)}\left(X_3\right)+f_{12}^{\left(2\right)}\left(X_1,X_2\right)+f_{13}^{\left(2\right)}\left(X_1,X_3\right)+f_{23}^{\left(2\right)}\left(X_2,X_3\right)+f_{123}^{\left(2\right)}\left(X_1,X_2,X_3\right)\\ & {-f_2^{\left(2\right)}\left(\alpha \right)-f_3^{\left(2\right)}\left(\beta \right)-f_{12}^{\left(2\right)}\left(X_1,\alpha \right)-f_{13}^{\left(2\right)}\left(X_1,\beta \right)-f_{23}^{\left(2\right)}\left(\alpha ,\beta \right)-f_{123}^{\left(2\right)}\left(X_1,\alpha ,\beta \right)\}}^2],\end{array}

where \mathrm{t}\mathrm{r}\left[{\left\{l_2\left(\mathit{Y},{\tilde{\mathit{Y}}}_1\right)\right\}}^2\right] is evaluated as follows

(51)\begin{array}{rl}\mathrm{t}\mathrm{r}\left[{\left\{l_2\left(\mathit{Y},{\tilde{\mathit{Y}}}_1\right)\right\}}^2\right]=& D_2+D_3+D_{12}+D_{13}+D_{23}+D_{123}+{\left\{f_2^{\left(1\right)}\left(\alpha \right)\right\}}^2+{\left\{f_2^{\left(2\right)}\left(\alpha \right)\right\}}^2+{\left\{f_3^{\left(1\right)}\left(\beta \right)\right\}}^2+{\left\{f_3^{\left(2\right)}\left(\beta \right)\right\}}^2\\ & +\mathrm{E}\left[{\left\{f_{12}^{\left(1\right)}\left(X_1,\alpha \right)\right\}}^2\right]+\mathrm{E}\left[{\left\{f_{12}^{\left(2\right)}\left(X_1,\alpha \right)\right\}}^2\right]+\mathrm{E}\left[{\left\{f_{13}^{\left(1\right)}\left(X_1,\beta \right)\right\}}^2\right]+\mathrm{E}\left[{\left\{f_{13}^{\left(2\right)}\left(X_1,\beta \right)\right\}}^2\right]\\ & +\mathrm{E}\left[{\left\{f_{23}^{\left(1\right)}\left(\alpha ,\beta \right)\right\}}^2\right]+\mathrm{E}\left[{\left\{f_{23}^{\left(2\right)}\left(\alpha ,\beta \right)\right\}}^2\right]+\mathrm{E}\left[{\left\{f_{123}^{\left(1\right)}\left(X_1,\alpha ,\beta \right)\right\}}^2\right]+\mathrm{E}\left[{\left\{f_{123}^{\left(2\right)}\left(X_1,\alpha ,\beta \right)\right\}}^2\right].\end{array}

Now to formulate global first order sensitivity index of input variable X_1, evaluate \mathrm{E}\left[\mathrm{t}\mathrm{r}\left[{\left\{l_2\left(\mathit{Y},{\tilde{\mathit{Y}}}_1\right)\right\}}^2\right]\right] by treating \alpha and \beta as random variables with PDF, p_{X_2}\left(x_2\right) and p_{X_3}\left(x_3\right), respectively. Therefore, Eq. (51) reduces to

(52)\mathrm{E}\left[\mathrm{t}\mathrm{r}\left[{\left\{l_2\left(\mathit{Y},{\tilde{\mathit{Y}}}_1\right)\right\}}^2\right]\right]=2\left(D_2+D_3+D_{12}+D_{13}+D_{23}+D_{123}\right).

In a similar manner, global main sensitivity indices with respect to input variables X_2 and X_3 can be evaluated based on l₂ norm. Using Eq. (39), the above expression can be written as follows

(53)\mathrm{E}\left[\mathrm{t}\mathrm{r}\left[{\left\{l_2\left(\mathit{Y},{\tilde{\mathit{Y}}}_1\right)\right\}}^2\right]\right]=2D^{\mathit{Y}}\left(1-S_1\right),

where Eq. (53) shows the relationship between the main sensitivity index based on l₂ norm and Sobol’s first order index with respect to the input variable X_1. In general, the relationship between main sensitivity index based on l₂ norm and Sobol’s first order sensitivity index is shown in Eq. (54).

(54)S_i=1-\frac{\mathrm{E}\left[\mathrm{t}\mathrm{r}\left[{\left\{l_2\left(\mathit{Y},{\tilde{\mathit{Y}}}_i\right)\right\}}^2\right]\right]}{2D^{\mathit{Y}}};i=1,2,3.

4.1 Remarks

1) The equivalence between multivariate Sobol’s and l₂ norm-based analyses can be shown only if the variables of the input vector \mathit{X} are independent.

2) The orthogonal properties of the summands of Sobol’s function as specified in Eqs. (6) to (8) hold true only if the variables of \mathit{X} are independent.

3) When the variables of \mathit{X} are dependent, the global total and main sensitivity indices based on l₂ norm can be evaluated with respect to the input random variable X_i using {\displaystyle \frac{\mathrm{E}\left[\mathrm{t}\mathrm{r}\left[{\left\{l_2\left(\mathit{Y},{\mathit{Y}}_i\right)\right\}}^2\right]\right]}{2D^{\mathit{Y}}}} and {\displaystyle \frac{\mathrm{E}\left[\mathrm{t}\mathrm{r}\left[{\left\{l_2\left(\mathit{Y},{\tilde{\mathit{Y}}}_i\right)\right\}}^2\right]\right]}{2D^{\mathit{Y}}}}, respectively. A higher value of total sensitivity index (or a lower value of main sensitivity index) indicates a more significant sensitivity impact or vice versa.

4) The sensitivity analyses of altered and reference models by fixing the input random variable(s) at nominal value(s) give the notion of global-local sensitivity indices. Using Eq. (44), the global-local total sensitivity index of X_1 can be evaluated by fixing it as \mu_. Similarly, the global-local main sensitivity index of X_1 can be evaluated by fixing X_2\mathrm{a}\mathrm{n}\mathrm{d}{X}_3 to \alpha \mathrm{a}\mathrm{n}\mathrm{d}\beta respectively as mentioned in Eq. (51).

5 Model distance-based sensitivity analysis for grouped variables

This section proposes performing multivariate global sensitivity analysis for grouped variables using model distance-based approach. Consider a reference model, \mathit{Y}=f\left(\mathit{X}\right). The variables in the input vector \mathit{X} are taken to be partitioned into two groups, {\mathit{X}}_{\mathit{u}} and {\mathit{X}}_{\sim \mathit{u}}. Here, \mathit{u} is a subset of \left\{1,2,\dots ,n\right\} and \sim \mathit{u} represents the complementary set of \mathit{u}. Consider another hypothetical model, where all the variables of {\mathit{X}}_{\sim \mathit{u}} are treated as uncertain. The variables in the group {\mathit{X}}_{\mathit{u}} are fixed at its nominal value. Let this model be called as altered model-1 and denoted by {\mathit{Y}}_{\mathit{u}}. The difference between the reference model, \mathit{Y} and the altered model-1, {\mathit{Y}}_{\mathit{u}} represents the total sensitivity index of the group {\mathit{X}}_{\mathit{u}}. This can be evaluated by comparing the joint PDFs of the models using probability distance measures such as Hellinger distance, Kullback–Leibler divergence, and l₂ norm.

Similarly, an altered model-2 can be defined where except for the group {\mathit{X}}_{\mathit{u}}, variables in the group {\mathit{X}}_{\sim \mathit{u}} are deterministically fixed at their respective nominal values. Let the response of this altered model-2 is denoted by {\tilde{\mathit{Y}}}_{\mathit{u}}. The difference between the reference model, \mathit{Y} and the altered model-2, {\tilde{\mathit{Y}}}_{\mathit{u}} represents the main sensitivity index of the group {\mathit{X}}_{\mathit{u}}.

5.1 Remarks

1) Following the notion discussed in Section 4, the global total sensitivity index of the grouped variable {\mathit{X}}_{\mathit{u}}, based on l_2 norm can be evaluated using {\displaystyle \frac{\mathrm{E}\left[\mathrm{t}\mathrm{r}\left[{\left\{l_2\left(\mathit{Y},{\mathit{Y}}_{\mathit{u}}\right)\right\}}^2\right]\right]}{2D^{\mathit{Y}}}}. The relationship between global total sensitivity index of {\mathit{X}}_{\mathit{u}}, based on l₂ norm and Sobol’s analysis is given by T_{{\mathit{X}}_{\mathit{u}}}={\displaystyle \frac{\mathrm{E}\left[\mathrm{t}\mathrm{r}\left[{\left\{l_2\left(\mathit{Y},{\mathit{Y}}_{\mathit{u}}\right)\right\}}^2\right]\right]}{2D^{\mathit{Y}}}}. This is valid when the variables of \mathit{X} are treated as independent. Here, T_{{\mathit{X}}_{\mathit{u}}}denotes Sobol’s total sensitivity index of {\mathit{X}}_{\mathit{u}} [56,57].

2) Similarly, the global main sensitivity index of the grouped variable {\mathit{X}}_{\mathit{u}}, based on l_2 norm can be evaluated using {\displaystyle \frac{\mathrm{E}\left[\mathrm{t}\mathrm{r}\left[{\left\{l_2\left(\mathit{Y},{\tilde{\mathit{Y}}}_{\mathit{u}}\right)\right\}}^2\right]\right]}{2D^{\mathit{Y}}}}. The relationship between global main sensitivity index of {\mathit{X}}_{\mathit{u}}, based on l₂ norm and Sobol’s analysis is given by S_{{\mathit{X}}_{\mathit{u}}}=1-{\displaystyle \frac{\mathrm{E}\left[\mathrm{t}\mathrm{r}\left[{\left\{l_2\left(\mathit{Y},{\tilde{\mathit{Y}}}_{\mathit{u}}\right)\right\}}^2\right]\right]}{2D^{\mathit{Y}}}}. This is valid when the variables of \mathit{X} are treated as independent. Here, S_{{\mathit{X}}_{\mathit{u}}} denotes Sobol’s first order sensitivity index of {\mathit{X}}_{\mathit{u}}.

3) If T_{{\mathit{X}}_{\mathit{u}}}=0, then it indicates that {\mathit{X}}_{\mathit{u}} group of input variables, does not contribute to the variability of the response \mathit{Y} either directly or by its interactions with other group of input variables.

4) If S_{{\mathit{X}}_{\mathit{u}}}=0, then it implies that {\mathit{X}}_{\mathit{u}} group of input variables, does not directly contributing to the variability of the response \mathit{Y}.

5) If T_{{\mathit{X}}_{\mathit{u}}}=S_{{\mathit{X}}_{\mathit{u}}}=1, then it implies that \mathit{Y} is solely dependent on {\mathit{X}}_{\mathit{u}} only.

6 Estimation of multivariate global sensitivity indices

This section discusses the procedure for estimating multivariate global total and main sensitivity indices using brute-force MCS. The procedure for computing the global-local sensitivity indices is discussed is steps I to IX, followed by the procedure for computing the global sensitivity indices.

Step I. Generate N samples of input random vector \mathit{X} from the joint PDF, p_{\mathit{X}}\left(x_1,x_2,\dots ,x_n\right) using MCS. Let these samples be denoted by {\left[x_1^{\left(j\right)},x_2^{\left(j\right)},\dots ,x_n^{\left(j\right)}\right]}^{\mathrm{t}};j=1,2,\dots ,N.

Step II. Evaluate the response of the reference model, \mathit{Y}={\left[Y^{\left(1\right)},Y^{\left(2\right)},\dots ,Y^{\left(s\right)}\right]}^{\mathrm{t}} using the samples simulated from step (I). Let the realizations of Y be denoted by {\left[y_j^{\left(1\right)},y_j^{\left(2\right)},\dots ,y_j^{\left(s\right)}\right]}^{\mathrm{t}}; j=1,2,\dots ,N. Here, y_j^{\left(s\right)} represents the j th realization of the random variable, Y^{\left(s\right)} belonging to the random vector, Y.

Step III. Evaluate the response of the altered model-1, {\mathit{Y}}_i={\left[Y_i^{\left(1\right)},Y_i^{\left(2\right)},\dots ,Y_i^{\left(s\right)}\right]}^{\mathrm{t}};i=1,2,\dots ,n by fixing X_i={\alpha }_i. Then, {\left[y_{i_j}^{\left(1\right)},y_{i_j}^{\left(2\right)},\dots ,y_{i_j}^{\left(s\right)}\right]}^{\mathrm{t}}=f(x_1^{(j)},x_2^{(j)},\dots ,{x_{i-1}^{(j)},{\alpha }_i,x}_{i+1}^{(j)},\dots ,x_n^{(j)});j=1,2,\dots ,N. Here, y_{i_j}^{\left(s\right)} represents the j th realization of the random variable, Y_i^{\left(s\right)} belonging to the random vector, {\mathit{Y}}_i.

Step IV. Evaluate the response of the altered model-2, {\tilde{\mathit{Y}}}_i={\left[{\tilde{Y}}_i^{\left(1\right)},{\tilde{Y}}_i^{\left(2\right)},\dots ,{\tilde{Y}}_i^{\left(s\right)}\right]}^{\mathrm{t}};i=1,2,\dots ,n by maintaining the uncertainty of X_{\mathrm{i}}, where all other input random variables are made deterministic by fixing them to constant values, {\alpha }_j;j=1,2,\dots ,n;j\ne i. Then, {\left[{\tilde{y}}_{i_j}^{\left(1\right)},{\tilde{y}}_{i_j}^{\left(2\right)},\dots ,{\tilde{y}}_{i_j}^{\left(s\right)}\right]}^{\mathrm{t}}=f({\alpha }_1,{\alpha }_2,\dots ,{\alpha }_{\mathrm{i}-1},x_i^{\left(j\right)},{\alpha }_{\mathrm{i}+1},\dots ,{\alpha }_n);j=1,2,\dots ,N. Here, {\tilde{y}}_{i_j}^{\left(s\right)} represents the jth realization of the random variable, {\tilde{Y}}_i^{\left(s\right)} belonging to the random vector, {\tilde{\mathit{Y}}}_i.

Step V. Estimate the joint PDF of the models \mathit{Y} and {\mathit{Y}}_i, denoted by {\hat{p}}_{\mathit{Y}}\left(y^{\left(1\right)},y^{\left(2\right)},\dots ,y^{\left(s\right)}\right) and {\hat{p}}_{{\mathit{Y}}_i}\left(y_i^{\left(1\right)},y_i^{\left(2\right)},\dots ,y_i^{\left(s\right)}\right) respectively, using a multivariate kernel density estimator (KDE) [58].

Step VI. The multivariate global-local total sensitivity index with respect to the input random variable, X_i based on Hellinger distance can be estimated by evaluating the distance between the joint PDFs of \mathit{Y} and {\mathit{Y}}_i, as shown in Eq. (57). Here, {\hat{p}}_{\mathit{Y}}\left(y_i^{\left(1\right)},y_i^{\left(2\right)},\dots ,y_i^{\left(s\right)}\right) and {\hat{p}}_{{\mathit{Y}}_i}\left(y^{\left(1\right)},y^{\left(2\right)},\dots ,y^{\left(s\right)}\right) are evaluated using an appropriate interpolation method.

(55)k_1=\underset{-\mathrm{\infty }}{\stackrel{+\mathrm{\infty }}{\iint }}\cdots \underset{-\mathrm{\infty }}{\stackrel{+\mathrm{\infty }}{\int }}\sqrt{\left({\hat{p}}_{\mathit{Y}}\left(y^{\left(1\right)},y^{\left(2\right)},\dots ,y^{\left(s\right)}\right){\hat{p}}_{{\mathit{Y}}_i}\left(y^{\left(1\right)},y^{\left(2\right)},\dots ,y^{\left(s\right)}\right)\right)}\prod _{r=1}^s\mathrm{d}{y}^{\left(r\right)},

(56)k_2=\underset{-\mathrm{\infty }}{\stackrel{+\mathrm{\infty }}{\iint }}\cdots \underset{-\mathrm{\infty }}{\stackrel{+\mathrm{\infty }}{\int }}\sqrt{\left({\hat{p}}_{\mathit{Y}}\left(y_i^{\left(1\right)},y_i^{\left(2\right)},\dots ,y_i^{\left(s\right)}\right){\hat{p}}_{{\mathit{Y}}_i}\left(y_i^{\left(1\right)},y_i^{\left(2\right)},\dots ,y_i^{\left(s\right)}\right)\right)}\prod _{r=1}^s\mathrm{d}{y}_i^{\left(r\right)},

(57){\delta }_{\mathrm{H}}\left({\alpha }_i\right)=1-0.5\left(k_1+k_2\right).

Step VII. The multivariate global-local total sensitivity index with respect to the input random variable, X_i based on Kullback–Leibler divergence can be estimated by evaluating the distance between the joint PDFs of \mathit{Y} and {\mathit{Y}}_i, as shown in Eq. (58).

(58)\begin{array}{rl}{\delta }_{\mathrm{K}}\left({\alpha }_i\right)=& \underset{-\mathrm{\infty }}{\stackrel{+\mathrm{\infty }}{\iint }}\cdots \underset{-\mathrm{\infty }}{\stackrel{+\mathrm{\infty }}{\int }}{\hat{p}}_{\mathit{Y}}\left(y^{\left(1\right)},y^{\left(2\right)},\dots ,y^{\left(s\right)}\right)\\ & \cdot \mathrm{l}\mathrm{o}\mathrm{g}\frac{{\hat{p}}_{\mathit{Y}}\left(y^{\left(1\right)},y^{\left(2\right)},\dots ,y^{\left(s\right)}\right)}{{\hat{p}}_{{\mathit{Y}}_i}\left(y^{\left(1\right)},y^{\left(2\right)},\dots ,y^{\left(s\right)}\right)}\prod _{r=1}^s\mathrm{d}{y}^{\left(r\right)}.\end{array}

Step VIII. The multivariate global-local total sensitivity index with respect to the input random variable, X_{\mathrm{i}} based on l₂ norm can be estimated as follows

(59){\hat{l}}_2\left(\mathit{Y},{\mathit{Y}}_i\right)={\displaystyle \frac{{\displaystyle \frac{1}{N}}{\displaystyle \sum _{r=1}^s{\displaystyle \sum _{j=1}^N{\left(y_j^{\left(r\right)}-y_{i_j}^{\left(r\right)}\right)}^2}}}{2D^{\mathit{Y}}}}.

Step IX. Repeat the steps (III) to (VIII) for all other input random variables, X_i;i=1,2,\dots ,n. In a similar manner, the multivariate global-local main sensitivity index with respect to the input random variable X_i can also be evaluated by estimating the distance between the joint PDFs of the reference model, \mathit{Y} and the altered model-2, {\tilde{\mathit{Y}}}_i.

The total and main sensitivity indices estimated in the previous steps (I) to (IX), depend upon the deterministic values, {\alpha }_i;i=1,2,\dots ,n and hence, these indices are considered as global-local sensitivity indices. Further, the global sensitivity indices can be estimated by repeatedly performing the steps I to IX for different values of {\alpha }_i;i=1,2,\dots ,n. The average of these estimated indices yields the global total and main sensitivity indices with respect to the input random variable X_i.

7 Illustration

To exemplify the proposed model distance-based multivariate sensitivity analysis, two engineering problems are examined in this illustration. The first example deals with the quantification of uncertainty in the deflections of a steel truss, and the second example entails quantifying the uncertainty in the displacements of a two-story steel frame subjected to ground motion.

7.1 Example 1-static analysis of a truss system

Consider a simply supported steel truss as shown in Fig.2. In this study, the structure is subjected to uniformly distributed loads acting in both horizontal and vertical directions, as denoted by q_{\mathrm{H}} and q_{\mathrm{V}}, respectively. The point loads, denoted by P_i;i=1,2,\dots ,5 are acting on this structure as shown in Fig.2.

Two responses considered in this example for performing the multivariate global sensitivity analysis and these are the vertical deflections at nodes, D and G, denoted by {\delta }_{\mathrm{D}} and {\delta }_{\mathrm{G}} respectively. The deflections are evaluated using Eq. (60). The axial rigidities of the inclined and horizontal members are denoted by {AE}_1 and {AE}_2, respectively.

\begin{array}{rl}{\delta }_{\mathrm{D}}=& \frac{1}{A{E}_1}(46.52q_{\mathrm{H}}+1.95q_{\mathrm{V}}-0.86P_1+2.58P_2\\ & +6.03P_3+0.86P_4+4.31P_5)+\frac{1}{A{E}_2}(18.21q_{\mathrm{H}}+0.65q_{\mathrm{V}}\\ & -1.58P_1+2.17P_2+2.45P_3+0.01P_4+3.46P_5),\end{array}

(60)\begin{array}{rl}{\delta }_{\mathrm{G}}=& \frac{1}{A{E}_1}(92.95q_{\mathrm{H}}+0.01q_{\mathrm{V}}+5.16P_1+7.75P_2\\ & +2.58P_3+10.33P_4+5.16P_5)+\frac{1}{A{E}_2}(51.99q_{\mathrm{H}}+2.59q_{\mathrm{V}}\\ & +0.58P_1+5.19P_2+2.89P_3+3.47P_4+5.19P_5).\end{array}

Two separate cases are considered in this example.

(1) Case 1: Here, all the input variables representing the load variables are modeled as independent random variables, and variables representing the axial rigidities are treated as deterministic; and

(2) Case 2: Here, all the input variables representing both axial rigidities and load variables are modeled as dependent random variables.

7.1.1 Case 1

For Case 1, the axial rigidities of the inclined members and horizontal members are taken as 4.73 × 10⁶ kN and 1.99 × 10⁶ kN, respectively. The statistical properties of the uncertain load variables are detailed in Tab.2. The uncertain load variables are assumed to be normally distributed with a standard deviation of 10%. The global sensitivity analysis is performed on the uncertain deflections using the model distance-based approach. The brute-force MCS is used to evaluate the sensitivity indices of the input variables, and the results are compared with the exact solutions obtained as shown in Tab.3. Here, the numbers enclosed in parentheses signify the ranks of the input variables. In this case, a sample size of 10³ is used for MCS.

The evaluation of exact indices is as follows. The total sensitivity index with respect to a specific input variable is evaluated by comparing the joint PDF for reference and altered models. The deflection expressions representing this reference and altered models using Eq. (60) simplify into a linear expression over the load variables, when both axial rigidities are modeled as deterministic variables. This linear nature of the response functions ensures that, when input variables being normally distributed, both reference and altered models emerge as normally distributed random vectors. Then, the sensitivity index can be evaluated by comparing the multivariate normal distributions of these random vectors. For this purpose, the closed-form expressions available for Hellinger distance [59] and Kullback–Leibler divergence between two multivariate normal distributions can be utilized. The evaluated sensitivity index is global-local in nature as it depends on the constant value of the input variable used in altered model. To formulate a global nature for the sensitivity index, perform an expectation operation over this global-local sensitivity index by treating the respective input variable as a random variable. This can be achieved by suitable numerical integration schemes. By following this approach, exact solutions based on Hellinger distance and Kullback–Leibler divergence are evaluated. The procedure for obtaining the exact solutions based on Sobol’s and l₂ norm-based analyses are detailed in supplementary materials.

From Tab.3, it can be observed that the rank ordering of the variables remains identical across various probability distance measures. The sensitivity measures presented in Tab.3 indicate that the variable q_{\mathrm{H}} has the most significant influence, and the variable q_{\mathrm{V}} has the least impact on the uncertainties in the deflection at nodes D and G. The results also imply that the exact solutions closely match with the results obtained from brute-force MCS across different probability distance measures. However, sampling fluctuations are noted in a few cases, which can be reduced by increasing the sample size for MCS. It is evident from the evaluated sensitivity indices of l₂ norm and Sobol’s analysis that they are identical. For a linear problem, when input variables are independent, the total and main sensitivity indices are exactly same for both l₂ norm-based analysis and Sobol’s analysis. The values tabulated in Tab.3 indicate that the total sensitivity indices estimated using Sobol’s analysis sum up to 0.9999, signifying that there is no functional interaction or dependency between the input variables. In the evaluation of main sensitivity indices based on Hellinger distance and Kullback–Leibler divergence using brute-force MCS, the correlation coefficient between the reference model and altered model-2 was observed to be above 0.95. This resulted in inconsistent values of the evaluated sensitivity indices which deem for requirement of efficient computational strategies.

7.1.2 Case 2

In the second case, the axial rigidities are modeled as random variables along with the uncertain load variables. The variables are also assumed to be correlated as specified by the correlation coefficient matrix given in Eq. (61).

(61)\rho =\left[\begin{array}{ccccccccc}1& 0.6& 0& 0& 0& 0& 0& 0& 0\\ 0.6& 1& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 1& 0.3& 0.7& 0.5& 0.3& 0& 0\\ 0& 0& 0.3& 1& 0.5& 0.6& 0.4& 0& 0\\ 0& 0& 0.7& 0.5& 1& 0.6& 0.6& 0& 0\\ 0& 0& 0.5& 0.6& 0.6& 1& 0.7& 0& 0\\ 0& 0& 0.3& 0.4& 0.6& 0.7& 1& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& 0.8\\ 0& 0& 0& 0& 0& 0& 0& 0.8& 1\end{array}\right].

The statistical properties of the uncertain axial rigidities are detailed in Tab.4. The axial rigidities are modeled as log-normally distributed random variables with a standard deviation of 10%.

The total and main sensitivity indices of the correlated input variables are evaluated using brute-force MCS and these results are tabulated in Tab.5 and Tab.6. Here, the numbers enclosed in parentheses signify the ranks of the input variables. In this case, a sample size of {10}^3 is used for MCS.

When the input variables are dependent, the global sensitivity indices based on l₂ norm can be evaluated using the expression discussed in clause 3 of Section 4.1. From Tab.5 and Tab.6, it can be observed that the ranking of the top four variables and the least four variables follow similar pattern across different probability distance measures. The evaluated sensitivity measures indicate that the variable {AE}_2 is found to be most significant variable toward the uncertainties in the deflection at nodes, D and G.

7.2 Example 2-dynamic analysis of a frame subjected to seismic excitation

Consider a two-story steel frame as shown in Fig.3, which is modeled as a multi-degree of freedom (MDOF) system. The mass of each floor is represented by m_i;i=1,2,3. The flexural rigidities of members AB, FE, and GH are denoted as {EI}_1 and for the members BC and ED, denoted by {EI}_2. The lengths of these members are denoted by L_1 and L_2, respectively. The damping ratio is assumed to be 0.5, denoted as \zeta. The frame is subjected to the north–south component of ground motion recorded in the 1940 El Centro earthquake, which had a magnitude of 7.1. This ground motion is denoted by x_{\mathrm{g}}\left(t\right) in Fig.3. For performing sensitivity analysis, the two responses considered in this example are maximum displacements observed at the first and second story of the steel frame.

The governing equation of motion for the MDOF system is expressed as follows

(62)\mathit{M}\ddot{\mathit{Y}}\left(t\right)+\mathit{C}\dot{\mathit{Y}}\left(t\right)+\mathit{K}\mathit{Y}\left(t\right)=\mathit{F}\left(t\right),

(63)\begin{array}{rl}& \mathit{M}=\left[\begin{array}{cc}{m}_1+m_2& 0\\ 0& m_3\end{array}\right],\mathit{C}=\alpha \mathit{M}+\beta \mathit{K},\mathit{K}=\left[\begin{array}{cc}{k}_1+k_2& -k_2\\ -k_2& k_2\end{array}\right],\\ & \mathit{Y}\left(t\right)=\left[\begin{array}{c}{y}_1\left(t\right)\\ y_2\left(t\right)\end{array}\right],\mathit{F}\left(t\right)=\left[\begin{array}{c}-\left(m_1+m_2\right){\ddot{x}}_{\mathrm{g}}\left(t\right)\\ -m_3{\ddot{x}}_{\mathrm{g}}\left(t\right)\end{array}\right],\end{array}

where \mathit{M},\mathit{C},\mathit{K},\mathit{Y}, and \mathit{F} represent the mass matrix, damping matrix, stiffness matrix, displacement vector, and force vector of the MDOF system, respectively. Let ‘\cdot’ denotes the time derivative, and k_1={\displaystyle \frac{36E{I}_1}{{L_1}^3}} and k_2={\displaystyle \frac{24E{I}_2}{{L_2}^3}} represent the equivalent stiffnesses of first and second stories of the frame, respectively. Let \alpha and \beta denote the Rayleigh damping coefficients, evaluated using {\displaystyle \frac{2\zeta {\omega }_1{\omega }_2}{{\omega }_1+{\omega }_2}} and {\displaystyle \frac{2\zeta }{{\omega }_1+{\omega }_2}}, respectively. Here, {\omega }_1 and {\omega }_2 represent the natural frequencies of the first and second stories of the frame, respectively. To determine the maximum story displacements, Eq. (62) is written in the state-space form, expressed as follows

(64)\dot{\mathit{z}}\left(t\right)=\mathit{A}\mathit{z}\left(t\right)+\mathit{B}\mathit{F}\left(t\right),

(65)\begin{array}{rl}\mathit{z}\left(t\right)=& {\left[\begin{array}{c}\mathit{Y}\left(t\right)\\ \dot{\mathit{Y}}\left(t\right)\end{array}\right]}_{4\times 1},\mathit{A}=\left[\begin{array}{cc}\mathbf{O}& \mathit{I}\\ -{\mathit{M}}^{-1}\mathit{K}& -{\mathit{M}}^{-1}\mathit{C}\end{array}\right],\\ & \mathrm{a}\mathrm{n}\mathrm{d}\mathit{B}=\left[\begin{array}{c}\mathbf{O}\\ {\mathit{M}}^{-1}\end{array}\right],\end{array}

where \mathit{z}\left(t\right) denotes the state vector, \mathbf{O} and \mathit{I} represent the null and identity matrices of size 2\times 2, respectively. The time history analysis, lasting 53.74 s, is performed with zero initial conditions (i.e., displacement and velocity), using ODE45 in MATLAB with time step of 0.02 s.

The input random vector is denoted by \mathit{X} (i.e., \mathit{X}={\left[m_1,m_2,m_3,{EI}_1,{EI}_2,L_1,L_2\right]}^{\mathrm{t}}). The statistical characteristics of the variables in this example are detailed in Tab.7. All the variables of \mathit{X} are assumed to be log-normally distributed with a standard deviation of 10%. While performing the multivariate sensitivity analysis for these displacements, two cases are considered: Case 1: all the variables of \mathit{X} are modeled as independent random variables; and Case 2: the random variables of \mathit{X} are grouped into three distinct groups. The components of \mathit{X} are partitioned into three groups denoted by {\mathit{X}}^{\left(i\right)};i=1,2,3 such that {\mathit{X}}^{\left(1\right)}={\left[m_1,m_2,m_3\right]}^{\mathrm{t}}, {\mathit{X}}^{\left(2\right)}=[{EI}_1,{EI}_2{]}^{\mathrm{t}}, and {\mathit{X}}^{\left(3\right)}={[L_1,L_2]}^{\mathrm{t}}.

The global sensitivity analysis is performed on the uncertain displacements using the model distance-based approach based on the probability distance measures such as l₂ norm, Hellinger distance, and Kullback–Leibler divergence. The brute-force MCS is used to evaluate the sensitivity indices of the ungrouped input variables, and the results are presented in Tab.8 and Tab.9. Here, the numbers enclosed in parentheses signify the ranks of the input variables. In this case, a sample size of {10}^3 is used for MCS.

It can be observed from Tab.8 and Tab.9 that the ranking of all variables remains similar across different probability distance measures. The sensitivity measures evaluated indicate that the variable L_2 is found to be the most significant variable and the variable m_1 is the least influential toward the uncertainties related to maximum story displacements.

In Case 2, the total and main sensitivity indices of the grouped input variables are evaluated using brute-force MCS and these results are tabulated in Tab.10 and Tab.11. Here, the numbers enclosed in parentheses signify the ranks of the input variables. A sample size of 2000 is used for MCS. The ranking of grouped variables remains consistent across different probability distance measures, as shown in Tab.10. Similarly, Tab.11 illustrates a consistent pattern in the ranking of the most significant group variables across these measures. The sensitivity measures evaluated indicate that {\mathit{X}}^{\left(3\right)} group of variables exerts most significant impact on the uncertainties associated with maximum story displacements.

8 Discussion

This section highlights the significant findings of the proposed model distance-based approach for performing multivariate global sensitivity analysis using probability distance measures.

1) Unlike the classical Sobol’s method, which performs sensitivity analysis through variance decomposition using second-order moments of the responses, the proposed framework utilizes joint PDFs instead. It incorporates probability distance measures such as Hellinger distance, Kullback–Leibler divergence, and l₂ norm for evaluating the global sensitivity indices.

2) Notably, the sensitivity indices based on l₂ norm exhibit equivalence to the classical multivariate Sobol’s sensitivity indices, as illustrated in Example 1 (Case 1). For instance, the sensitivity index for X_1 is 0.7023 in both Sobol’s and l₂ norm-based analyses. A similar pattern is observed for the indices evaluated for the remaining input variables (Tab.3).

3) The proposed framework is capable of quantifying uncertainties in engineering systems with correlated input variables that follow non-Gaussian distributions, whereas the classical Sobol’s analysis is limited to handling only independent inputs.

4) Additionally, the proposed framework can perform multivariate global sensitivity analysis for a system with grouped variables as illustrated in Example 2 (Case 2).

5) The rank ordering of the input variables based on the (total and main) sensitivity indices, remains consistent across different probability distance measures. This has been observed in the case of grouped variables also.

6) In this study, the brute-force MCS approach used to estimate the sensitivity indices, demonstrated consistency with the exact solutions. For instance, in Example 1 (Case 1), the exact sensitivity index for variable X_2, computed using l₂ norm, is 0.0027, closely matching with the value obtained via MCS for a sample size of 10³. Similarly, the exact sensitivity indices based on Hellinger distance and Kullback–Leibler divergence are 0.0007 and 0.0027, respectively, both consistent with the MCS results (Tab.3).

However, sampling fluctuations are observed in a few cases. This issue can be reduced by increasing the sample size during MCS. For instance, in Example 1 (Case 1), the exact sensitivity index for the variable X_1 computed using l₂ norm is 0.7023. With a sample size of 10³, MCS estimate is 0.6660 (as reported in Tab.3). But when the sample size is increased to 10⁴, the estimated value converges to 0.7032, closely matching with the exact solution.

7) A major drawback of using the brute-force MCS within the proposed framework is the significant computational burden it entails. For an engineering system with n input variables, calculating the global total and main sensitivity indices requires \left(N+2{nN}^2\right) calls on the response model, where N represents the number of Monte Carlo samples. This indicates that the total number of simulations needed to compute both the total and main sensitivity indices increases substantially with larger sample sizes and/(or) with a greater number of input variables. This challenge is particularly pronounced in complex engineering problems that necessitates extensive computational effort to evaluate the response of model.

8) The limitation of multivariate KDE in high-dimensional output spaces is known, as its accuracy decreases due to the curse of dimensionality. To address this, a copula-based approach [60] can be used, which offers a scalable alternative that improves the accuracy of joint PDF estimation for high-dimensional problems.

9 Conclusions

This paper defines a novel concept of model distance-based multivariate sensitivity analysis, using probability distance measures such as l₂ norm, Hellinger distance, and Kullback–Leibler divergence, to quantify the uncertainties associated with multiple responses in an engineering system with respect to the uncertainties in the system parameters. The proposed methodology performs sensitivity analysis by using joint PDFs of the responses, instead of relying on their second-order moments as done in Sobol’s analysis. While Sobol’s analysis is effective for independent random variables, the proposed analysis is versatile and capable of handling situations where input random variables are dependent and non-Gaussian. This paper also demonstrates the equivalence between the l₂ norm-based analysis and Sobol’s multivariate global sensitivity analysis. Therefore, the multivariate global sensitivity analysis based on l₂ norm can be considered as a generalization of Sobol’s analysis.

Illustrative examples showcase the effectiveness of the proposed methodology in finding the most and least significant input random variables, which contributes to the uncertainties in system responses. The sensitivity measures evaluated through brute force MCS in Example 1 (Case 1) closely matches the exact solutions. This ensures the accuracy of brute force MCS approach across different probability distance measures. Moreover, the proposed methodology can effectively handle both correlated and non-Gaussian random variables. This has been illustrated in Example 1 (Case 2). The proposed methodology can also perform multivariate sensitivity analysis for systems with grouped variables, further enhancing the overall capability of the framework as illustrated in Example 2. In these illustrations, the consistent rank ordering of input variables, based on total and main sensitivity indices across different probability distance measures, demonstrates the effectiveness of the proposed methodology. Hence, the proposed methodology of model distance-based multivariate global sensitivity analysis can be applied to engineering systems. For instance, in structural health monitoring of large infrastructure such as bridges, multivariate global sensitivity analysis enables the evaluation of sensitivity indices for various input factors such as boundary conditions, material properties, load conditions, and environmental factors, considering their influence on multiple structural responses, including displacements, stress distributions, and strain measurements. The proposed framework quantifies the influence of input factors on the critical structural responses, addressing inherent uncertainties in engineering systems and facilitating the identification of both the most and least significant input variables.

However, the complexity of the MCS approach increases with high-dimensional problems and with larger sample sizes. This presents a significant computational challenge. To overcome this limitation, future work could explore new methodologies aimed at developing effective alternatives to brute-force MCS techniques through variance reduction schemes in complex engineering systems. Additionally, improving sampling techniques, such as Latin Hypercube Sampling [61], Quasi-Monte Carlo methods [62], or stratified sampling, could reduce the number of simulations required while maintaining accuracy. These ideas would help in enhancing computational efficiency and extend the applicability of the proposed framework to more intricate engineering problems.

References

[1] Der KiureghianADitlevsenO. Aleatory or epistemic? Does it matter? Structural Safety, 2009, 31(2): 105–112

[2] Chutia R , Mahanta S , Datta D . Uncertainty modelling of atmospheric dispersion by stochastic response surface method under aleatory and epistemic uncertainties. Sadhana, 2014, 39(2): 467–485

[3] Uusitalo L , Lehikoinen A , Helle I , Myrberg K . An overview of methods to evaluate uncertainty of deterministic models in decision support. Environmental Modelling and Software, 2015, 63: 24–31

[4] SaltelliARattoMAndresTCampolongoFCariboniJGatelliDTarantolaS. Global Sensitivity Analysis: The Primer. New Jersey: John Wiley and Sons, 2008

[5] Song X , Sun G , Li Q . Sensitivity analysis and reliability-based design optimization for high-strength steel tailor welded thin-walled structures under crashworthiness. Thin-walled Structures, 2016, 109: 132–142

[6] Ionescu-Bujor M , Cacuci D G . A comparative review of sensitivity and uncertainty analysis of large-scale systems: Deterministic methods. Nuclear Science and Engineering, 2004, 147(3): 189–203

[7] SmithR C. Uncertainty Quantification: Theory, Implementation, and Applications. Philadelphia, PA: Society for Industrial and Applied Mathematics, 2013

[8] Cotter S C . A screening design for factorial experiments with interactions. Biometrika, 1979, 66(2): 317–320

[9] Morris M D . Factorial sampling plans for preliminary computational experiments. Technometrics, 1991, 33(2): 161–174

[10] Sobol I M . Sensitivity estimates for nonlinear mathematical models. Mathematical Modelling and Computational Experiments, 1993, 4: 407–414

[11] Vu-Bac N , Lahmer T , Keitel H , Zhao J , Zhuang X , Rabczuk T . Stochastic predictions of bulk properties of amorphous polyethylene based on molecular dynamics simulations. Mechanics of Materials, 2014, 68: 70–84

[12] Vu-Bac N , Lahmer T , Zhang Y , Zhuang X , Rabczuk T . Stochastic predictions of interfacial characteristic of polymeric nanocomposites (PNCs). Composites. Part B, Engineering, 2014, 59: 80–95

[13] Vu-Bac N , Silani M , Lahmer T , Zhuang X , Rabczuk T . A unified framework for stochastic predictions of mechanical properties of polymeric nanocomposites. Computational Materials Science, 2015, 96: 520–535

[14] Borgonovo E . A new uncertainty importance measure. Reliability Engineering and System Safety, 2007, 92(6): 771–784

[15] Sobol I M . Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates. Mathematics and Computers in Simulation, 2001, 55(1–3): 271–280

[16] Vuillod B , Montemurro M , Panettieri E , Hallo L . A comparison between Sobol’s indices and Shapley’s effect for global sensitivity analysis of systems with independent input variables. Reliability Engineering and System Safety, 2023, 234: 109177

[17] Vu-Bac N , Rafiee R , Zhuang X , Lahmer T , Rabczuk T . Uncertainty quantification for multiscale modeling of polymer nanocomposites with correlated parameters. Composites. Part B, Engineering, 2015, 68: 446–464

[18] Vu-Bac N , Lahmer T , Zhuang X , Nguyen-Thoi T , Rabczuk T . A software framework for probabilistic sensitivity analysis for computationally expensive models. Advances in Engineering Software, 2016, 100: 19–31

[19] Vu-Bac N , Nguyen A H , Luong V H . Uncertainty analysis of static fatigue of Hi-Nicalon bundles. Engineering Analysis with Boundary Elements, 2024, 167: 105862

[20] Borgonovo E , Plischke E . Sensitivity analysis: A review of recent advances. European Journal of Operational Research, 2016, 248(3): 869–887

[21] Saltelli A , Jakeman A , Razavi S , Wu Q . Sensitivity analysis: A discipline coming of age. Environmental Modelling and Software, 2021, 146: 105226

[22] GrigoriuM. Stochastic Calculus: Applications in Science and Engineering. Boston: Birkhauser, 2002

[23] LutesL DSarkaniS. Random Vibration: Analysis of Structural and Mechanical Systems. Boston: Elsevier, 2004

[24] MelchersR EBeckA T. Structural Reliability Analysis and Prediction. New Jersey: John Wiley and Sons, 2018

[25] Kucherenko S , Rodriguez-Fernandez M , Pantelides C , Shah N . Monte Carlo evaluation of derivative-based global sensitivity measures. Reliability Engineering and System Safety, 2009, 94(7): 1135–1148

[26] Borgonovo E , Tarantola S . Moment independent and variance-based sensitivity analysis with correlations: An application to the stability of a chemical reactor. International Journal of Chemical Kinetics, 2008, 40(11): 687–698

[27] Marino S , Hogue I B , Ray C J , Kirschner D E . A methodology for performing global uncertainty and sensitivity analysis in systems biology. Journal of Theoretical Biology, 2008, 254(1): 178–196

[28] Khorashadi Zadeh F , Nossent J , Sarrazin F , Pianosi F , van Griensven A , Wagener T , Bauwens W . Comparison of variance-based and moment-independent global sensitivity analysis approaches by application to the SWAT model. Environmental Modelling and Software, 2017, 91: 210–222

[29] Pianosi F , Wagener T . A simple and efficient method for global sensitivity analysis based on cumulative distribution functions. Environmental Modelling and Software, 2015, 67: 1–11

[30] Chun M H , Han S J , Tak S J . An uncertainty importance measure using a distance metric for the change in a cumulative distribution function. Reliability Engineering and System Safety, 2000, 70(3): 313–321

[31] Greegar G , Manohar C S . Global response sensitivity analysis using probability distance measures and generalization of Sobol’s analysis. Probabilistic Engineering Mechanics, 2015, 41: 21–33

[32] Greegar G , Manohar C S . Global response sensitivity analysis of uncertain structures. Structural Safety, 2016, 58: 94–104

[33] Nandi S , Singh T . Global sensitivity analysis measures based on statistical distances. International Journal for Uncertainty Quantification, 2021, 11(6): 1–30

[34] Abhinav S , Manohar C S . Global response sensitivity analysis of randomly excited dynamic structures. Journal of Engineering Mechanics, 2016, 142(3): 04015094

[35] Greegar G , Abhinav S , Manohar C S . Model distance–based global–local response-sensitivity indexes for randomly inhomogeneous structures under stochastic excitations. ASCE-ASME Journal of Risk and Uncertainty in Engineering Systems. Part A, Civil Engineering, 2018, 4(3): 05018002

[36] Li X , Lei Y , Liu L . Efficient moment-independent sensitivity analysis of uncertainties in seismic demand of bridges based on a novel four-point-estimate method. Applied Sciences, 2021, 11(21): 10405

[37] Campbell K , McKay M D , Williams B J . Sensitivity analysis when model outputs are functions. Reliability Engineering and System Safety, 2006, 91(10–11): 1468–1472

[38] Zhang K , Lu Z , Cheng K , Wang L , Guo Y . Global sensitivity analysis for multivariate output model and dynamic models. Reliability Engineering and System Safety, 2020, 204: 107195

[39] Anderson M J . A new method for non-parametric multivariate analysis of variance. Austral Ecology, 2001, 26(1): 32–46

[40] Gamboa F , Janon A , Klein T , Lagnoux A . Sensitivity indices for multivariate outputs. Comptes Rendus Mathematique, 2013, 351(7–8): 307–310

[41] Lamboni M , Monod H , Makowski D . Multivariate sensitivity analysis to measure global contribution of input factors in dynamic models. Reliability Engineering and System Safety, 2011, 96(4): 450–459

[42] Xiao S , Lu Z , Xu L . Multivariate sensitivity analysis based on the direction of eigen space through principal component analysis. Reliability Engineering and System Safety, 2017, 165: 1–10

[43] Xiao S , Lu Z , Wang P . Multivariate global sensitivity analysis based on distance components decomposition. Risk Analysis, 2018, 38(12): 2703–2721

[44] Xiao S , Lu Z , Wang P . Global sensitivity analysis based on distance correlation for structural systems with multivariate output. Engineering Structures, 2018, 167: 74–83

[45] Li L , Lu Z , Wu D . A new kind of sensitivity index for multivariate output. Reliability Engineering and System Safety, 2016, 147: 123–131

[46] Rizzo M L , Szekely G . J. Energy distance. Wiley interdisciplinary reviews: Computational statistics, 2016, 8(1): 27–38

[47] Xiao S , Lu Z , Wang P . Multivariate global sensitivity analysis for dynamic models based on energy distance. Structural and Multidisciplinary Optimization, 2018, 57(1): 279–291

[48] Liu F , Wei P , Tang C , Wang P , Yue Z . Global sensitivity analysis for multivariate outputs based on multiple response Gaussian process model. Reliability Engineering and System Safety, 2019, 189: 287–298

[49] Sun X , Choi Y Y , Choi J I . Global sensitivity analysis for multivariate outputs using polynomial chaos-based surrogate models. Applied Mathematical Modelling, 2020, 82: 867–887

[50] Liu Y , Li L , Chang Z , Wang P . Kriging-based analytical technique for global sensitivity analysis of systems with multivariate output. Structural and Multidisciplinary Optimization, 2022, 65(10): 288

[51] Chen L , Huang H . Global sensitivity analysis for multivariate outputs using generalized RBF-PCE metamodel enhanced by variance-based sequential sampling. Applied Mathematical Modelling, 2024, 126: 381–404

[52] TangWJiangWYangXZhengQ. Global sensitivity analysis of high-speed train dynamics system with high-dimensional inputs and multiple outputs. Vehicle System Dynamics, 2025: 1–21

[53] ChianiLBorgonovoEPlischkeETavoniM. Global sensitivity analysis of integrated assessment models with multivariate outputs. Risk Analysis, 2025: 1

[54] Mazo G , Tournier L . An inference method for global sensitivity analysis. Technometrics, 2025, 67(2): 270–282

[55] RachevS TKlebanovL BStoyanovS VFabozziF. The Methods of Distances in the Theory of Probability and Statistics. New York: Springer, 2013

[56] Gamboa F , Janon A , Klein T , Lagnoux A . Sensitivity analysis for multidimensional and functional outputs. Electronic Journal of Statistics, 2014, 8(1): 575–603

[57] Broto B , Bachoc F , Depecker M , Martinez J M . Sensitivity indices for independent groups of variables. Mathematics and Computers in Simulation, 2019, 163: 19–31

[58] Botev Z I , Grotowski J F , Kroese D P . Kernel density estimation via diffusion. The Annals of Statistics, 2010, 38(5): 2916–2957

[59] PardoL. Statistical Inference Based on Divergence Measures. Boca Raton, FL: Chapman and Hall/CRC, 2018

[60] NelsenR B. An Introduction to Copulas. New York: Springer, 2006

[61] Vu-Bac N , Zhuang X , Rabczuk T . Uncertainty quantification for mechanical properties of polyethylene based on fully atomistic model. Materials, 2019, 12(21): 3613

[62] Nguyen T C , Vu-Bac N , Budarapu P R . A Kriging-based uncertainty quantification for fracture assessment. International Journal of Computational Methods, 2025, 22(1): 2450045

RIGHTS & PERMISSIONS

Higher Education Press