PCESensitivityAnalysis¶

Purpose¶

The goal of the driver is to perform a variance based sensitivity analysis of a scalar function $f: \mathcal{X} \rightarrow \mathbb{R}$ or a vectorial function $\mathbf{f}: \mathcal{X} \rightarrow \mathbb{R}^m$ with respect to variations of a parameter vector $\mathbf{p}\in\mathcal{X}\subset\mathbb{R}^d$ . To this end, the function is decomposed into poylnomials using a polynomial chaos expansion (PCE).

Depending on the defined random distribution of each parameter (see parameter distribution), PCE uses different polynomials to obtain a good $L^2$ convergence.

Distribution of random variable	Polynomials $\{\phi_j\}_{j\in\mathbb{N}_0}$
Uniform	Legendre
Gauss	Hermite
Gamma	Laguerre
Beta	Jacobi

Using the default setting, the number of of expansion terms is chosen automatically based on a variance analysis of the model vector $\mathbf{y}(\mathbf{p})$ .

The PCE can be used to predict the function and its derivatives with respect to each entry of the parameter vector $\mathbf{p}$ . Predictions can be computed by calling study.predict().

There are two typical application scenarios of a variance based sensitivity analysis:

The scalar function $f(\mathbf{p})$ measures the performance of a device. The sensitivity analysis helps to quantify the average performance for random perturbations of the parameter vector $\mathbf{p}$ . The analysis also quantifies which parameter combinations influence the performance .
The vectorial function $\mathbf{f}(\mathbf{p})$ describes a list of physical observables for parameter reconstruction. The analysis allows to identify vector entries with high sensitivity for specific parameters.

The parameter vector $\mathbf{p}$ consists of $d$ entries, where each entry $k$ has a specific random distribution $\rho_k(p_k)$ (see parameter distribution). Hence, the scalcar function $f(\mathbf{p})$ is also a random variable (vectorial outputs are analysed by multiple independent sensitivity analyses).

The variance of this random variable can analysed in terms of Sobol coefficients $S_{i_1,\dots,i_s}$ . Each Sobol coefficient determines the fraction of the total variation that stems from the joint random variation of the parameters $p_{i_1},\dots,p_{i_s}$ . The sum of all Sobol coefficients is 1,

$\sum_{i=1}^d S_i + \sum_{i<j}^s S_{i j} + \dots + S_{1 2 \dots d} = 1$

For example, if $S_{1 3} = 0.1$ , then 10% of the variance of $f(\mathbf{p})$ stem from the variance of an additive part $f_{1 3}(p_1,p_3)$ of the function that only depends on the random variations of parameter $p_1$ and $p_3$ . For more details on the definition of Sobol coefficients, see the section Additional Information.

Usage Example¶

addpath(fullfile(getenv('JCMROOT'), 'ThirdPartySupport', 'Matlab'));
client = jcmwave_optimizer_client();

% Definition of the search domain
domain = {
    struct('name','x1', 'type','continuous', 'domain',[-5,5]),...
    struct('name','x2', 'type','continuous', 'domain',[-5,5]),...
    struct('name','x3', 'type','continuous', 'domain',[-5,5]),...
};

% Creation of the study object with study_id 'example'
study = client.create_study('domain',domain,  ...
                'driver','PCESensitivityAnalysis',...
                'name','PCESensitivityAnalysis example', ...
                'study_id','PCESensitivityAnalysis_example');

% Definition of the objective function (Ishigami function)
a = 7; b = 0.1; %parameters of the Ishigami function
function obs = objective(sample)
  x1 = sample.x1;
  x2 = sample.x2;
  x3 = sample.x3;

  obs = study.new_observation();
  %objective
  obs = obs.add(sin(x1)+ a*sin(x2)^2 + b*x3^4*sin(x1));

end

% Define random distribution of parameters
distribution = {
    struct('name','x1', 'distribution','uniform', 'domain',[-pi,pi]),...
    struct('name','x2', 'distribution','uniform', 'domain',[-pi,pi]),...
    struct('name','x3', 'distribution','uniform', 'domain',[-pi,pi]) ...
};
% Set study parameters
study.set_parameters('max_iter',150, 'distribution',distribution);

% Run the study
while(not(study.is_done))
    sug = study.get_suggestion();
    obs = objective(sug.sample);
    study.add_observation(obs, sug.id);
end

%Analytic mean, variance and Sobol coefficients for the Ishigami function
mean = a/2;
variance = 1/2 + a^2/8 + b*pi^4/5 + b^2*pi^8/18;
sobol_values = struct(...
    'x1', 0.5*(1+b*pi^4/5)^2/variance,...
    'x2', (a^2/8)/variance,...
    'x1_x3', b^2*pi^8/2*(1/9-1/25)/variance...
);

driver_info = study.driver_info();
fprintf('\nMean %.3f Exact %.3f', driver_info.mean, mean);
fprintf('\nVariance %.3f Exact %.3f', driver_info.variance, variance);
fns = fieldnames(driver_info.sobol_indices);
for i=1:length(fns)
    sobol_index = driver_info.sobol_indices.(fns{i});
    val_driver = driver_info.sobol_coefficients.(fns{i});
    try
        val_exact = sobol_values.(fns{i});
    catch
        val_exact = 0;
    end
    fprintf('\nSobol index %s: %.3f Exact %.3f',mat2str(sobol_index),...
                                                val_driver,val_exact);
end

Parameters¶

The following parameters can be set by calling, e.g.

study.set_parameters('example_parameter1',[1,2,3], 'example_parameter2',true);

max_iter (int):	Maximum number of evaluations of the objective function (default: inf)

max_time (int):	Maximum run time in seconds (default: inf)

num_parallel (int):
	Number of parallel observations of the objective function (default: 1)

distribution (list):

Definition of random distribution for each parameter in the format of a list. All continuous parameters with unspecified distribution are assumed to be uniformely distributed in the parameter domain. Fixed and discrete parameters are not random parameters. The value of discrete parameters defaults to the first listed value. (default: None)

Example:

{struct("name","param1", "distribution","normal", "mean",1.0, "variance",2.0),
 struct("name","param2", "distribution","uniform", "domain",[-5,5]),
 struct("name","param3", "distribution","gamma", "alpha",1.2, "beta",0.5),
 struct("name","param4", "distribution","fixed", "value",7.0),
 struct("name","param5", "distribution","beta", "alpha",1.5, "beta",0.8)}

sampling_strategy (string):
	Sampling strategy of parameter values. standard: Sampler for cartesian product samples of given parameters under algebraic constraints. WLS: Sampler for weighted Least-Squares sampling as described by Cohen and Migliorati. (default: WLS) (options: [‘standard’, ‘WLS’])

polynomial_degree (int):
	Maximum degree of the polynomial chaos expansion. Can be also a list of degrees for every parameter. If not set, the polynomial degree is iteratively adapted. (default: None)

limit_polynomial_degree (int):
	Maximum sum of polynomial degree of the chaos expansion. (default: 2147483647)

optimization_step (int):
	Number of iterations before the polynomial degrees of the expansion are optimized. Only applies if the polynomial degree is determined automatically (polynomial_degree=None). (default: 10)

optimization_step_max (int):
	Maximum number of iterations after which the polynomial degrees of the expansion are optimized. Only applies if the polynomial degree is determined automatically (polynomial_degree=None). (default: 1000)

pce_update_step (int):
	Number of iterations before the PCE is updated from the samples drawn so far. At each update, e.g. the values of the Sobol coefficients are updated. (default: 5)

num_test_samples (int):
	Number of samples that are drawn in order to estimate the prediction error of the PCE. (default: 10)

min_prediction_error (float):
	Stopping criterium. Minimum error of the prediction of the polynomial chaos expansion for the last 5 aqcuisitions. (default: 1e-06)

target_cond (float):
	Target condition number of the information matrix. Only applies if the polynomial degree is determined automatically (polynomial_degree=None). (default: 500)

Driver-specific Information¶

A struct with following parameters can be retrieved by calling study.driver_info().

sobol_indices:	Dictionary mapping the names of Sobol indices (e.g. x1_x3) to the list of the corresponding parameter indices (e.g. [1,3]).

sobol_coefficients:
	Dictionary mapping the name of the Sobol indices (e.g. x1_x3) to the value of the Sobol coefficients.

input_covariance:
	Covariance matrix between entries of the vectorial input function $\mathbf{f}(\mathbf{p})$ , i.e. $C_ij = \mathbb{E}\left[(f_i(\mathbf{p})-\overline{f}_i)(f_j(\mathbf{p})-\overline{f}_j)\right]$

predicition_error:
	Maximal predicion error of test data.

num_expansion_terms:
	Number of PCE expansion terms.

mean:	Mean of the objective functions under parameter uncertainties.

variance:	Variance of the objective functions under parameter uncertainties.

Additional Information¶

Sobol coefficients¶

Each function $f(\mathbf{p})$ can be decomposed into a unique sum called a Sobol decomposition,

$f(\mathbf{p}) = f_0 + \sum_{i=1}^d f_i(p_i) + \sum_{i<j}^{d} f_{ij}(p_i,p_j) + \cdots + f_{1,2,\dots,d}(p_1,p_2,\dots,p_d).$

Here, $f_0$ is a constant and each univariate integral over each term vanishes, i.e.

$\int f_{i_1,i_2,\dots,i_s}(p_{i_1},p_{i_2},\dots,p_{i_s}) \rho_k(p_k) \text{d}p_k = 0 \text{ for all } k=i_1,i_2,\dots,i_s.$

That is, all the terms in the functional decomposition are orthogonal. The variance of the function $f(\mathbf{p})$ is given as

$\text{Var}[f] = \int f^2(\mathbf{p}) \rho_1(p_1)\cdots\rho_d(p_d)\text{d}p_1\cdots\text{d}p_k - f_0^2.$

Due to the orthogonality, the variance can be decomposed into a sum of variances stemming from each term in the Sobol decomposition

$\text{Var}[f] = \sum_i \text{Var}[f_i] + \sum_{i<j} \text{Var}[f_ij] + \cdots + \text{Var}[f_{1 2\dots d}]$

where

$\text{Var}[f_{i_1,i_2,\dots,i_s}] = \int f_{i_1,i_2,\dots,i_s}^2(p_{i_1},p_{i_2},\dots,p_{i_s}) \rho_{i_1}(p_{i_1})\cdots\rho_{i_s}(p_{i_s})\text{d}p_{i_1}\cdots\text{d}p_{i_s}.$

The Sobol coefficient $S_{i_1,\dots,i_s}$ is defined as the fraction of the total variance that stems from the additive term $f_{i_1,i_2,\dots,i_s}$ , i.e.

$S_{i_1,\dots,i_s} = \frac{\text{Var}[f_{i_1,i_2,\dots,i_s}]}{\text{Var}[f]}.$