Solved – global sensitivity analysis input and output matrix

I am performing nonlinear analysis in ABAQUS software, with 5 input variables. For each realization of the input vector, I run the code and I get a scalar response value. I am using the Monte Carlo method to generate a sample for the input vector. I generated a Monte Carlo sample of size $N=500000$ for the input vector, and correspondigly I obtained a sample of the same size for the response by running the code. How can I perform a global sensitivity analysis (GSA) without a predefined function for it? Thanks

Consider a function of $p$ random variables

$$Y=f(X_1,dots,X_p)$$

When the output is neither a linear nor a monotone function of the outputs, but $X_1,dots,X_p$ are independent, and f is square-integrable, then a good way to estimate the global sensitivity of $Y$ to $X_1,dots,X_p$ is to use the Sobol indices. The definition of Sobol indices is related to the ANOVA decomposition(note this is a permalink). Assuming that you are only interested in the main effect indices and total order indices, and assuming that you are really going to run $N=5cdot10^6$ nonlinear FEM analyses and storing all the corresponding results (?), you can use fast99 from the R package sensitivity:

N <- 10^6 p <- 5 library(sensitivity) x <- fast99(model = NULL, factors = p, n = N, q = "qnorm", q.arg = list(mean = 0, sd = 1)) 

x is an object which, among the other things, stores the design you will need to run, i.e., the $M=pN$ realizations of the input vector, corresponding to the $M$ runs of your code. You can access them with

x$X

you can consult the docs of sensitivity for additional information. Once you have the $M$ runs of your code, and the corresponding vector of responses y, you can update the object x with

tell(x, y) 

Now, S <- x$D1/x$V and T <- 1- x$Dt/x$V are respectively the first-order indices (or main effect indices) and the total order indices.

fast99 is reasonably fast, as far as non-model based methods for the estimate of Sobol indices go. However, it has two drawbacks:

  1. it only gives you a point estimate, but not a confidence interval
  2. it can be inaccurate when the Sobol indices are small in magnitude w.r.t. 1, even though with millions of runs this won't be a problem.

An alternative to fast99, which is more accurate in the estimation of the smaller indices, but require thrice the number of runs, is sobolowen. This returns a Monte-Carlo based estimator of the Sobol indices, thus you get a confidence interval along with the point estimates. In this case you generate 3 design matrices this way:

X1 <- data.frame(matrix(runif(p * N), nrow = N)) X2 <- data.frame(matrix(runif(p * N), nrow = N)) X3 <- data.frame(matrix(runif(p * N), nrow = N)) 

and then, for 95%-confidence interval Monte Carlo estimates of the main effect and total order indices, you call

x <- sobolowen(model = NULL, X1, X2, X3, nboot = 100) 

As before, x$X contains your DOE. Update x with tell, and now x$S and x$T store respectively the first-order indices and the total order indices. Consult the vignette of sensitivity for many other estimators of the Sobol indices, with varying computational efforts.

PS note that here I assumed all your variables to be standard normal. I guess $X_1,dots,X_p$ are design parameters, whose distributions you assumed normal, but which are likely not standard normal. This means that to obtain your actual DOE, you need to transform each column of x$X by multiplying it for the standard deviation of variable $X_i$, and then add the mean of $X_i$.

EDIT: I just got the abstract of this paper in my feeds;

http://www.ans.org/pubs/journals/nt/a_38989

It discusses exactly your problem (though I assume their ABAQUS model will be more complex than yours), so you should give it a look. If you cannot access the paper from the above link, these two pdfs seem quite similar, even if a little older (1 year):

https://inldigitallibrary.inl.gov/sti/6468109.pdf

https://inldigitallibrary.inl.gov/sti/6799579.pdf

Similar Posts:

Rate this post

Leave a Comment