Polynomial Chaos Surrogate
This example explains how to train and utilize a PolynomialChaos surrogate model. For detailed information on the polynomial chaos method including the mathematics and implementation see PolynomialChaos and QuadratureSampler. For general information on training and evaluating a surrogate model see Training a surrogate model and Evaluating a surrogate model.
Overview
Polynomial chaos (PC) is a stochastic element method. Its formulation is similar to finite element method, but instead of discretizing spatial dimensions, PC discretizes random variables, i.e. uncertain parameters. Within the training phase, PC creates a functional of a quantity of interest (QoI) that is dependent on the random variables. This function consists of a finite sum of orthogonal polynomials with coefficients, these coefficients are what the training phase is attempting to compute. Another difference from finite elements is that these polynomials do not necessarily have compact support and increasing the terms in the sum increases the polynomial order. Please see PolynomialChaos for more details of the polynomial chaos method.
Model Problem
This example uses a one-dimensional heat conduction problem as the full-order model which has certain uncertain parameters. The model equation is as follows:
The quantities of interest are the average and maximum temperature:
Parameter Uncertainty
For demonstration, each of these parameters will have two types of probability distributions: Uniform () and Normal (). Where and are the max and minimum bounds of the uniform distribution, respectively. And and are the mean and standard deviation of the normal distribution, respectively.
The uncertain parameters for this model problem are:
| Parameter | Symbol | Uniform | Normal |
|---|---|---|---|
| Conductivity | |||
| Volumetric Heat Source | |||
| Domain Size | |||
| Right Boundary Temperature |
Analytical Solutions
This simple model problem has analytical descriptions for the field temperature, average temperature, and maximum temperature:
With the quadratic feature of the field temperature, using quadratic elements in the discretization will actually yield the exact solution.
Input File
Below is the input file used to solve the one-dimensional heat conduction model.
[Mesh]
type = GeneratedMesh
dim = 1
nx = 100
xmax = 1
elem_type = EDGE3
[]
[Variables]
[T]
order = SECOND
family = LAGRANGE
[]
[]
[Kernels]
[diffusion]
type = MatDiffusion
variable = T
diffusivity = k
[]
[source]
type = BodyForce
variable = T
value = 1.0
[]
[]
[Materials]
[conductivity]
type = GenericConstantMaterial
prop_names = k
prop_values = 2.0
[]
[]
[BCs]
[right]
type = DirichletBC
variable = T
boundary = right
value = 300
[]
[]
[Executioner]
type = Steady
solve_type = PJFNK
petsc_options_iname = '-pc_type -pc_hypre_type'
petsc_options_value = 'hypre boomeramg'
[]
[Postprocessors]
[avg]
type = AverageNodalVariableValue
variable = T
[]
[max]
type = NodalExtremeValue
variable = T
value_type = max
[]
[]
[Outputs]
[]
With this input the uncertain parameters are defined as:
Materials/conductivity/prop_valuesKernels/source/valueMesh/xmaxBCs/right/value
These values in the sub.i file are arbitrary since the stochastic master app will be modifying them.
Training
This section describes how to set up an input for a PolynomialChaosTrainer. There are three aspects that need to be known before setting up the input file.
A sub app needs to be built, defining the physics of the problem, this example uses a heat conduction model described in the previous section.
Uncertain parameters need to be identified with a defined probability distribution.
Quantities of interest need to be defined, for now, these come in the form of postprocessors from the sub app.
Omitting Solve
Any input file in MOOSE needs to include a Mesh, Variables, and Executioner block. However, the stochastic master app does not actually create or solve a system. So the StochasticToolsAction builds a minimal model to satisfy these requirements:
[StochasticTools]
[]
Parameter Uncertainty Distributions
A required aspect of the polynomial chaos method is a definition of the probability distribution of the uncertain parameters. These distributions define the type of polynomial used for the PC expansion. This example shows the two types of distributions available for PC training:
Uniform distribution for each parameter
[Distributions]
[k_dist]
type = Uniform
lower_bound = 1
upper_bound = 10
[]
[q_dist]
type = Uniform
lower_bound = 9000
upper_bound = 11000
[]
[L_dist]
type = Uniform
lower_bound = 0.01
upper_bound = 0.05
[]
[Tinf_dist]
type = Uniform
lower_bound = 290
upper_bound = 310
[]
[]
Normal distribution for each parameter
[Distributions]
[k_dist]
type = Normal
mean = 5
standard_deviation = 2
[]
[q_dist]
type = Normal
mean = 10000
standard_deviation = 500
[]
[L_dist]
type = Normal
mean = 0.03
standard_deviation = 0.01
[]
[Tinf_dist]
type = Normal
mean = 300
standard_deviation = 10
[]
[]
Training Sampler
Non-intrusive polynomial chaos training relies on sampling the full-order model (i.e. sub app) at a number of sample points. There are two general types of sampling techniques for PC: Monte Carlo and numerical quadrature.
Monte Carlo Sampling
Monte Carlo sampling is the most intuitive type of sampling and can be done using the MonteCarlo or LatinHypercube samplers in the stochastic tools module.
[Samplers]
[sample]
type = MonteCarlo
num_rows = 10000
distributions = 'k_dist q_dist L_dist Tinf_dist'
execute_on = PRE_MULTIAPP_SETUP
[]
[]
Quadrature Sampling
Because PC expansion uses very specific types of polynomials, using numerical quadrature is typically a more precise approach to training a PC model. The QuadratureSampler implements numerical quadrature for uniform and normal distributions. There are three different techniques for building a multidimensional quadrature grid: Cartesian grid, Smolyak sparse grid, and Clenshaw-Curtis sparse grid. These grids are specified by the "sparse_grid" parameter: none, smolyak, and clenshaw-curtis. Using sparse grids can significantly reduce the number of training points while retaining accuracy, see QuadratureSampler for more details. The number of sample points is ultimately defined by the "order" parameter. It is generally advisable that this parameter be set to the same value as the "order" parameter in PolynomialChaosTrainer.
[Samplers]
[sample]
type = Quadrature
order = 10
distributions = 'k_dist q_dist L_dist Tinf_dist'
execute_on = PRE_MULTIAPP_SETUP
[]
[]
Running Sub Application
Running the sub app and transferring data back and forth for PC is exactly the same as any other type of surrogate. See Training a surrogate model on more details on setting up this part of the input file.
[MultiApps]
[sub]
type = SamplerFullSolveMultiApp
input_files = sub.i
sampler = sample
[]
[]
[Controls]
[cmdline]
type = MultiAppSamplerControl
multi_app = sub
sampler = sample
param_names = 'Materials/conductivity/prop_values Kernels/source/value Mesh/xmax BCs/right/value'
[]
[]
[Transfers]
[data]
type = SamplerReporterTransfer
from_multi_app = sub
sampler = sample
stochastic_reporter = results
from_reporter = 'avg/value max/value'
[]
[]
[Reporters]
[results]
type = StochasticReporter
[]
[]
Trainers
The PC training object is defined within the Trainers block. The required parameters for a PC trainer are:
"distributions" specify the type of polynomials to use for the expansion, it is very import that these distributions match the distributions given to the sampler.
"sampler" is the object that will provide sample points which are given to the sub app during execution.
"response" specifies the result vector for storing the computed values.
"order" defines the maximum order of the PC expansion, this parameter ultimately defines the accuracy and complexity of the surrogate model.
[Trainers]
[poly_chaos_avg]
type = PolynomialChaosTrainer
execute_on = timestep_end
order = 10
regression_type = integration
distributions = 'k_dist q_dist L_dist Tinf_dist'
sampler = sample
response = results/data:avg:value
[]
[poly_chaos_max]
type = PolynomialChaosTrainer
execute_on = timestep_end
order = 10
regression_type = integration
distributions = 'k_dist q_dist L_dist Tinf_dist'
sampler = sample
response = results/data:max:value
[]
[]
Outputting Training Data
Outputting the data after training a PC trainer is exactly the same as outputting data for any other type of trainer. See Training a surrogate model on more details on setting up this part of the input file.
[Outputs]
file_base = poly_chaos_training
[out]
type = SurrogateTrainerOutput
trainers = 'poly_chaos_avg poly_chaos_max'
execute_on = FINAL
[]
[]
Example Input Files
Evaluation and Statistics
This section will go over how to set up an input file to evaluate and perform statistical analysis on a trained polynomial chaos surrogate. The polynomial chaos surrogate is unique in that it is able to compute statistical moments and sensitivities analytically. Specific postprocessors are available to compute these quantities without the need for sampling.
Omitting Solve
Any input file in MOOSE needs to include a Mesh, Variables, and Executioner block. However, the stochastic master app does not actually create or solve a system. So the StochasticToolsAction builds a minimal model to satisfy these requirements:
[StochasticTools]
[]
Evaluation
Evaluating a polynomial chaos surrogate model is exactly the same as any other type of surrogate, see Evaluating a surrogate model for more details.
Loading training data
[Surrogates]
[poly_chaos_avg]
type = PolynomialChaos
filename = 'poly_chaos_training_poly_chaos_avg.rd'
[]
[poly_chaos_max]
type = PolynomialChaos
filename = 'poly_chaos_training_poly_chaos_max.rd'
[]
[]
Defining sampler for evaluation – Uniform distributions
[Distributions]
[k_dist]
type = Uniform
lower_bound = 1
upper_bound = 10
[]
[q_dist]
type = Uniform
lower_bound = 9000
upper_bound = 11000
[]
[L_dist]
type = Uniform
lower_bound = 0.01
upper_bound = 0.05
[]
[Tinf_dist]
type = Uniform
lower_bound = 290
upper_bound = 310
[]
[]
[Samplers]
[sample]
type = MonteCarlo
num_rows = 100000
distributions = 'k_dist q_dist L_dist Tinf_dist'
execute_on = initial
[]
[]
Defining sampler for evaluation – Normal distributions
[Distributions]
[k_dist]
type = Normal
mean = 5
standard_deviation = 2
[]
[q_dist]
type = Normal
mean = 10000
standard_deviation = 500
[]
[L_dist]
type = Normal
mean = 0.03
standard_deviation = 0.01
[]
[Tinf_dist]
type = Normal
mean = 300
standard_deviation = 10
[]
[]
[Samplers]
[sample]
type = MonteCarlo
num_rows = 100000
distributions = 'k_dist q_dist L_dist Tinf_dist'
execute_on = initial
[]
[]
Evaluating surrogate with EvaluateSurrogate
[Reporters]
[samp]
type = EvaluateSurrogate
model = 'poly_chaos_avg poly_chaos_max'
sampler = sample
parallel_type = ROOT
[]
[]
Statistics and Sensitivities
PolynomialChaos has the ability to compute statistical moments, local sensitivity, and sobol indices analytically. These are computed using PolynomialChaosReporter. The implemented moments include mean, standard deviation, skewness, and kurtosis. The type of moment is specified by the "statistics" parameter, which can be a combination of mean, stddev, skewness, and kurtosis. This example computes mean and standard deviation. It should be noted that computing skewness and kurtosis can be very computationally demanding using the analytical technique. Therefore, it might be better to sample the surrogate model and compute these quantities with the result.
The points to compute the sensitivity can be defined explicitly by specifying the "local_sensitivity_points" parameter, and/or using a sampler by specifying the "local_sensitivity_sampler". It is vital that the columns of the points match with the distributions defined from the original trainer.
Sobol statistics are a metric of the global sensitivity of each parameter, or a combination of parameters, see SobolReporter for more details. These are computed by setting the "include_sobol" parameter to true. This will compute total, first-, and second-order Sobol indices.
[Reporters]
[stats]
type = PolynomialChaosReporter
pc_name = 'poly_chaos_avg poly_chaos_max'
statistics = 'mean stddev'
local_sensitivity_points = '5 10000 0.03 300; 5 10000 0.03 300'
include_sobol = true
[]
[]
Example Input Files
Results and Analysis
In this section, the results of training and evaluation of the surrogate model as described in the previous sections is shown. A convergence study of the polynomial chaos training, including analysis of sampling and polynomial order, is also provided.
Evaluation Results
Figure 1 and Figure 2 show the resulting probability distributions from sampling each surrogate with 100,000 points.
Figure 1: Temperature distributions with uniform parameter distribution
Figure 2: Temperature distribution with normal parameter distribution
Local Sensitivity Results
Table 1 shows the results of computing the local sensitivity of each parameter for every surrogate. For reference, the sensitivity is defined as
where is either or , is , , , or , and the point tested () is
Table 1: Local sensitivity results
| QoI | Distribution | Sensitivity | Sensitivity | Sensitivity | Sensitivity |
|---|---|---|---|---|---|
| Uniform | -0.001944 | 0.001985 | 0.003976 | 0.9977 | |
| Normal | -0.002009 | 0.001667 | 0.003303 | 0.9983 | |
| Uniform | -0.002918 | 0.002979 | 0.005966 | 0.9967 | |
| Normal | -0.002897 | 0.002736 | 0.005437 | 0.9973 |
Sobol Statistics Results
Figure 3 shows the results of computing Sobol statistics for average temperature and uniform parameter distributions (the other surrogates produced very similar statistics). In the plot, the x axis represents the first index subscript and the y axis represents the second. As can be seen, the statistics are symmetric, i.e. it does not matter the order of the index subscripts.
Figure 3: Average temperature Sobol statistics results uniform parameter distribution
Sampling Analysis
Here the effect of various sampling techniques on the training of a polynomial chaos surrogate is shown. Three different types of samplers were used: Monte Carlo, tensor quadrature, and Smolyak sparse quadrature. Monte Carlo was run with increasing number of samples and the quadrature was run with increasing order. Note, the polynomial expansion order was set constant. The performance metric is the error in the mean and standard deviation versus the number of total training points. Figure 4-Figure 7 show the results of this study for each quantity of interest and each parameter distribution. The uniform distributions (Figure 4 and Figure 5) have an analytical mean and standard deviation, which is why there is nice convergence. However, the normal distributions (Figure 6 and Figure 7) use numerical integration to compute the reference mean and standard deviation, which is why there is a more spurious convergence at low errors.
Figure 4: Sample convergence for moments of average temperature with uniform parameter distributions
Figure 5: Sample convergence for moments of maximum temperature with uniform parameter distributions
Figure 6: Sample convergence for moments of average temperature with normal parameter distributions
Figure 7: Sample convergence for moments of maximum temperature with normal parameter distributions
Polynomial Order Analysis
The effect of polynomial expansion order on the accuracy of the surrogate model is shown here. Again, the performance metric is the error in the mean and standard deviation versus the number of total training points. The maximum polynomial order is ranged from 3 to 9, where the tensor quadrature sampler order matches (which is generally advisable in any case). Figure 8 shows the result of the convergence study with a uniform distribution. Normal distribution was omitted since there is no analytical reference.
Figure 8: Polynomial order convergence for moments of maximum and average temperature with uniform parameter distributions