PCE METAMODELING: MULTIPLE OUTPUTS¶

This example showcases an application of polynomial chaos expansion (PCE) to a multi output model.

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 (7d8bd5ddc7274cf989d7b1dbba1cc7bd) started.

 uqpylab.sessions :: INFO     :: Reset successful.

Set the random seed for reproducibility¶

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

Computational model¶

Create a MODEL object from the function:

In [5]:
ModelOpts = {
    'Type': 'Model',
    'ModelFun': 'simply_supported_beam_9points.model',
    'isVectorized': 'true'
}

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', # beam width
        'Type': 'Lognormal',
        'Moments': [0.15, 0.0075] # (m)
        },
        {
        'Name': 'h', # beam height
        'Type': 'Lognormal',
        'Moments': [0.3, 0.015] # (m)
        },
        {
        'Name': 'L', # beam length
        'Type': 'Lognormal',
        'Moments': [5, 0.05] # (m)
        },
        {
        'Name': 'E', # Young's modulus
        'Type': 'Lognormal',
        'Moments': [3e10, 4.5e9] # (Pa)
        },
        {
        'Name': 'p', # uniform load
        'Type': 'Lognormal',
        'Moments': [1e4, 1e3] # (N/m)
        }]
}

myInput = uq.createInput(InputOpts)

Polynomial chaos expansion (PCE) metamodels¶

Select PCE as the metamodeling tool:

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

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

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

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

In [9]:
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.

Least-square methods rely on the evaluation of the model response on an experimental design. The following options configure UQLab to generate an experimental design of size 150 based on a latin hypercube sampling of the input model (also available: 'MC', 'Sobol', 'Halton'):

In [10]:
MetaOpts['ExpDesign'] = {
    "NSamples" : 120,
    "Sampling" : "LHS"
}

Create the LARS-based PCE metamodel:

In [11]:
myPCE = uq.createModel(MetaOpts)

 uqpylab.sessions :: INFO     :: Received intermediate compute request, function: simply_supported_beam_9points.model.

 uqpylab.sessions :: INFO     :: Carrying out local computation...

 uqpylab.sessions :: INFO     :: Local computation complete.

 uqpylab.sessions :: INFO     :: Starting transmission of intermediate compute results ((120, 9))...

 uqpylab.sessions :: INFO     :: Intermediate compute results sent.

Processing . done!

Validation of the results¶

The deflections $V(s_i)$ at the nine points are plotted for three realizations of the random inputs. Relative length units are used for comparison, because $L$ is one of the random inputs.

Create and evaluate the validation set¶

Generate validation samples:

In [12]:
Nval = 3
Xval = uq.getSample(myInput, Nval)

Evaluate the original "simply supported beam" model at the validation set points:

In [13]:
Yval = uq.evalModel(myModel, Xval)

Evaluate the PCE metamodels at the same points:

In [14]:
YPC = uq.evalModel(myPCE, Xval)

Create plots¶

For each sample points of the validation set $\mathbf{x}^{(i)}$, the simply supported beam deflection $\mathcal{M}(\mathbf{x}^{(i)})$ is compared against the one predicted by the PCE metamodel $\mathrm{\mathcal{M}^{PCE}(\mathbf{x}^{(i)})}$:

In [15]:
li = np.linspace(0,1,11,endpoint=True)
legend_text = []

for idx in np.arange(Nval):
    x=li
    y=np.concatenate([[0],YPC[idx,:],[0]])
    pp = plt.plot(x,y,':',linewidth=1.0, c=uq_colors[idx])
    legend_text.append('$\\mathrm{\\mathcal{M}}(\\mathbf{x}^{('+str(idx+1)+')}$)')

    y=np.concatenate([[0],Yval[idx,:],[0]])
    plt.plot(x,y,'-',linewidth=1.0,c=uq_colors[idx])
    # Add text for the plot legend
    legend_text.append('$\\mathrm{\\mathcal{M}^{PCE}}(\\mathbf{x}^{('+str(idx+1)+')}$)')

plt.xlabel('$\\mathrm{L}_{rel}$')
plt.ylabel('V(m)')
plt.grid(True)
plt.legend(legend_text)
plt.show()
No description has been provided for this image

Terminate the remote UQCloud session¶

In [16]:
mySession.quit()

 uqpylab.sessions :: INFO     :: Session 7d8bd5ddc7274cf989d7b1dbba1cc7bd terminated.
Out[16]:
True