PCE METAMODELING: ESTIMATION OF STATISTICAL MOMENTS¶

This example showcases an application of polynomial chaos expansion (PCE) to the estimation of the statistical moments of a computational model with random inputs. The model is a simply supported beam model that computes deflections at mid-span of a beam subjected to a uniform load. The convergence behavior of the PCE-based and pure Monte-Carlo-based estimation of the statistical moments is compared.

Package imports¶

In [1]:
from uqpylab import sessions, display_util
import numpy as np
import matplotlib.pyplot as plt

Initialize common plotting parameters¶

In [2]:
display_util.load_plt_defaults()
uq_colors = display_util.get_uq_color_order()

Start a remote UQCloud session¶

In [3]:
# Start the session
mySession = sessions.cloud()
# (Optional) Get a convenient handle to the command line interface
uq = mySession.cli
# Reset the session
mySession.reset()

Processing .
.
 done!

 uqpylab.sessions :: INFO     :: This is UQ[py]Lab, version 1.00, running on https://uqcloud.ethz.ch. 
                                 UQ[py]Lab is free software, published under the open source BSD 3-clause license.
                                 To request special permissions, please contact:
                                  - Stefano Marelli (marelli@ibk.baug.ethz.ch).
                                 A new session (3ff59606a67741dc879027ba85da660b) started.

 uqpylab.sessions :: INFO     :: Reset successful.

Set the random seed for reproducibility¶

In [4]:
uq.rng(100,'twister');

Computational model¶

The simply supported beam model is shown in the following figure:

SimplySupportedBeam

The forward model computes the deflection of the beam $V$ at mid-span location, which reads:

$$ V = \frac{5pL^4}{32Ebh^3} $$

This computation is carried out by the function simply_supported_beam(X) supplied with UQCloud. The input variables of this function are gathered into the $N \times M$ matrix X, where $N$ and $M$ are the numbers of the realizations and input variables, respectively. The input variables are given in the following order:

  1. $b$: beam width ($\text{m}$)
  2. $h$: beam height ($\text{m}$)
  3. $L$: beam length ($\text{m}$)
  4. $E$: Young's modulus ($\text{MPa}$)
  5. $p$: uniform load ($\text{kN/m}$)

Create a MODEL object from the function:

In [5]:
ModelOpts = {
    'Type': 'Model',
    'ModelFun': 'simply_supported_beam.model'
}
myModel = uq.createModel(ModelOpts)

Probabilistic input model¶

The simply supported beam has five inputs, modeled by independent lognormal random variables. The detailed model is given in the following table:

| Variable | Description | Distribution | Mean | Std.Deviation | |:--------:|:------------|:------------:|:------:|:--------------: | $b$ | Beam width | Lognormal | $0.15$ $\text{m}$ | $7.5$ $\text{mm}$ | | $h$ | Beam height | Lognormal | $0.3$ $\text{m}$ | $15$ $\text{mm}$| | $L$ | Length | Lognormal | $5$ $\text{m}$ | $50$ $\text{mm}$ | | $E$ | Young's modulus | Lognormal | $30000$ $\text{MPa}$ | $4500$ $\text{MPa}$ | | $p$ | Uniform load | Lognormal | $10$ $\text{kN/m}$ | $2$ $\text{kN/m}$ |

Define an INPUT object with the following marginals:

In [6]:
InputOpts = {
    "Marginals": [
        {
            "Name": "b",
            "Type": "Lognormal",
            "Moments": [0.15, 0.0075]
        },
        {
            "Name": "h",
            "Type": "Lognormal",
            "Moments": [0.3, 0.015]
        },
        {
            "Name": "L",
            "Type": "Lognormal",
            "Moments": [5, 0.05]
        },
        {
            "Name": "E",
            "Type": "Lognormal",
            "Moments": [3e10, 4.5e9]
        },
        {
            "Name": "p",
            "Type": "Lognormal",
            "Moments": [1e4, 2e3]
        }
    ]
}

Create the INPUT object:

In [7]:
myInput = uq.createInput(InputOpts)

Polynomial chaos expansion (PCE) metamodel¶

Select PCE as the metamodeling tool:

In [8]:
MetaOpts = {
    'Type': 'Metamodel',
    'MetaType': 'PCE'
}

Select the sparse-favoring least-square minimization LARS for the coefficient calculation strategy:

In [9]:
MetaOpts['Method'] = 'LARS'

Specify a sparse truncation scheme (hyperbolic norm with $q = 0.75$):

In [10]:
MetaOpts['TruncOptions'] = {'qNorm': 0.75}

Specify the range of the degrees to be compared by the adaptive algorithm:

In [11]:
MetaOpts['Degree'] = np.arange(2,11).tolist()

The degree with the lowest Leave-One-Out cross-validation (LOO) error estimator is chosen as the final metamodel.

Monte-Carlo- vs. PCE-based estimation of moments¶

The moments of a PCE can be analytically calculated from its coefficients without additional sampling. In this section, the mean and standard deviation estimates obtained from Monte-Carlo (MC) samples of increasing size are compared to the ones based on the PCE of the same samples.

The simply supported beam model with lognormal input distributions allows for the explicit calculation of mean and standard deviation, which are provided for reference:

In [12]:
mean_exact = 0.008368320689566
std_exact = 0.002538676671701

Specify a set of experimental design sizes to test:

In [13]:
NED = [50, 100, 150, 200, 500]

Initialize the arrays in which the results will be stored:

In [14]:
mean_MC = np.zeros(len(NED))
std_MC = np.zeros(len(NED))
mean_PCE = np.zeros(len(NED))
std_PCE = np.zeros(len(NED))

Loop over the experimental design sizes. In each step, generate an experimental design, create a PCE metamodel, and compute the corresponding mean and standard deviation:

In [15]:
for i, n in enumerate(NED):
    # Get a sample from the probabilistic input model with LHS sampling
    X_ED = uq.getSample(myInput, n, 'LHS')

    # Evaluate the full model on the current experimental design
    Y_ED = uq.evalModel(myModel, X_ED)
    
    # Calculate the moments of the experimental design
    mean_MC[i] = np.mean(Y_ED)
    std_MC[i] = np.std(Y_ED)
    
    # Use the sample as the experimental design for the PCE
    MetaOpts['ExpDesign'] = {
        'X': X_ED.tolist(), 
        'Y': Y_ED.tolist()
    }
    
    # Create and add the PCE model to UQCloud
    myPCE = uq.createModel(MetaOpts)
    
    # Calculate the mean and standard deviation from the PCE coefficients
    mean_PCE[i] = myPCE['PCE']['Moments']['Mean']
    std_PCE[i] = np.sqrt(myPCE['PCE']['Moments']['Var'])

Convergence plots¶

Plot the convergence of the mean estimates from MCS and PCE:

In [16]:
plt.semilogy(NED, np.abs(mean_MC/mean_exact - 1), linestyle="-", color=[0, 0.2980, 0.6000])
plt.semilogy(NED, np.abs(std_MC/std_exact - 1), linestyle="--", color=[0, 0.2980, 0.6000])
plt.semilogy(NED, np.abs(mean_PCE/mean_exact - 1), linestyle="-", color=[0.9412, 0.5294, 0])
plt.semilogy(NED, np.abs(std_PCE/std_exact - 1), linestyle="--", color=[0.9412, 0.5294, 0])
plt.ylim([1e-8,1])
plt.xlabel('N')
legend_text = ["$\\hat{\\mu}^{MC}_Y/\\mu_Y - 1$",
               "$\\hat{\\sigma}^{MC}_Y/\\mu_Y - 1$",
               "$\\hat{\\mu}^{PCE}_Y/\\sigma_Y - 1$",
               "$\\hat{\\sigma}^{PCE}_Y/\\sigma_Y - 1$"]
plt.legend(legend_text, frameon=False, loc="best")
plt.grid(True, which="both")
plt.show()
No description has been provided for this image

Terminate the remote UQCloud session¶

In [17]:
mySession.quit()

 uqpylab.sessions :: INFO     :: Session 3ff59606a67741dc879027ba85da660b terminated.
Out[17]:
True