KRIGING METAMODELING: MULTIPLE OUTPUTS

This example showcases an application of Kriging to the metamodeling of a simply supported beam model with multiple outputs. The model computes the deflections at several points along the length of the beam subjected to a uniform random load.

Package imports

from uqpylab import sessions, display_util
import numpy as np 
import matplotlib.pyplot as plt

# Initialize common plotting parameters
display_util.load_plt_defaults()
uq_colors = display_util.get_uq_color_order()

Start a remote UQCloud session

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

 uqpylab.sessions :: INFO     :: A new session (cc4bf4eea91946e79f6b7b1920cbeb2b) started.
 uqpylab.sessions :: INFO     :: Reset successful.

Set the random seed for reproducibility

uq.rng(50,'twister');

Computational model

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

The (negative) deflection of the beam at any longitudinal coordinate $s$ is given by:

$$V(s) = -\frac{p \,s (L^3 - 2\, s^2 L + s^3) }{2E b h^3}$$

This computation is carried out by the function simply_supported_beam_9points. The function evaluates the inputs gathered in the $N \times M$ matrix X, where $N$ and $M$ are the numbers of realizations and inputs, respectively. The inputs are given in the following order:

• $b$: beam width $(m)$
• $h$: beam height $(m)$
• $L$: beam length $(m)$
• $E$: Young's modulus $(Pa)$
• $p$: uniform load $(N/m)$

The function returns the beam deflection $V(s_i)$ at nine equally-spaced points along the length $s_i = i \cdot L/10, \; i=1,\ldots,9.$

Create a MODEL object from the simply_supported_beam_9points function:

ModelOpts = {
    'Type': 'Model',
    'ModelFun': 'simply_supported_beam_9points.model',
    'isVectorized': 'true'
}

myModel = uq.createModel(ModelOpts)

Probabilistic input model

The simply supported beam model has five independent input parameters modeled by lognormal random variables. The parameters of the distributions are given in the following table:

Variable	Description	Distribution	Mean	Std. deviation
b	Beam width	Lognormal	0.15 m	7.5 mm
h	Beam height	Lognormal	0.3 m	15 mm
L	Length	Lognormal	5 m	50 mm
E	Young modulus	Lognormal	30000 MPa	4500 MPa
p	Uniform load	Lognormal	10 kN/m	2 kN/m

Define an INPUT object with the following marginals:

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)

Kriging metamodel

Select the Kriging metamodeling tool:

MetaOpts = {
    "Type": "Metamodel",
    "MetaType": "Kriging"
}

Assign the previously created INPUT and MODEL objects to the metamodel options:

MetaOpts["Input"] = myInput["Name"]
MetaOpts["FullModel"] = myModel["Name"]

If an INPUT object and a MODEL object are specified, an experimental design is automatically generated and the corresponding model responses are computed.

Specify the sampling strategy and the number of sample points for the experimental design:

MetaOpts["ExpDesign"] = {
    "Sampling" : "LHS",
    "NSamples" : 150
}

Set a linear trend:

MetaOpts["Trend"] = {
    "Type" : "linear"
}

Use the exponential correlation family:

MetaOpts["Corr"] = {
    "Family" : "exponential"
}

Use maximum-likelihood-based hyperparameter estimation:

MetaOpts["EstimMethod"] = "ML"

Use the hybrid genetic algorithm (GA) method to solve the hyperparameters estimation problem and specify the number of generations and population size:

MetaOpts["Optim"] = {
    "Method" : "HGA",
    "MaxIter" : 50,
    "HGA" : {
        "nPop" : 50
    }
}

Create the Kriging metamodel:

myKriging = 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 ((150, 9))...
 uqpylab.sessions :: INFO     :: Intermediate compute results sent.

Validation

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  is one of the random inputs.

Create and evaluate a validation sample

Generate a validation sample:

Nval = 3
Xval = uq.getSample(myInput, Nval)

# Evaluate the full computational model at the validation sample:
Yval = uq.evalModel(myModel, Xval)

# Evaluate the Kriging metamodel at the same sample points:
YKRG = uq.evalModel(myKriging, Xval)

Visualize results

For each sample point 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 Kriging metamodel $\mathcal{M}^{KRG}(\mathbf{x}^{(i)})$ (its mean):

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

for idx in np.arange(0,Nval):
    x=li
    y=np.concatenate([[0],YKRG[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}^{KRG}}(\mathbf{x}^{('+str(idx+1)+')}$)')

plt.xlabel('$\mathrm{L}_{rel}$')
plt.ylabel('V(m)')
plt.grid(True)
plt.legend(legend_text);
plt.show()

Terminate the remote UQCloud session

mySession.quit()