Note

Go to the end to download the full example code.

# Create a polynomial chaos metamodel by integration on the cantilever beam¶

In this example, we create a polynomial chaos metamodel by integration on the cantilever beam example. We choose to evaluate the coefficients of the chaos decomposition by integration using various kinds of design of experiment:

Gauss product

Latin hypercube sampling

Quasi Monte Carlo with a Sobol’ sequence

We will compare the results obtained on each design.

```
from openturns.usecases import cantilever_beam
import openturns as ot
import openturns.viewer as otv
ot.Log.Show(ot.Log.NONE)
```

We first load the model from the usecases module :

```
cb = cantilever_beam.CantileverBeam()
```

In this example we consider all marginals independent.
They are defined in the `CantileverBeam`

class:

```
dist_E = cb.E
dist_F = cb.F
dist_L = cb.L
dist_I = cb.II
distribution = cb.independentDistribution
```

We load the model.

```
dim_input = cb.dim # dimension of the input
dim_output = 1 # dimension of the output
g = cb.model
```

## Create a polynomial chaos decomposition¶

We create the multivariate polynomial basis by tensorization of the univariate polynomials and the default linear enumerate rule.

```
multivariateBasis = ot.OrthogonalProductPolynomialFactory(
[dist_E, dist_F, dist_L, dist_I]
)
```

In this case, we select P using the
`getBasisSizeFromTotalDegree()`

method,
so that all polynomials with total degree lower or equal to 5 are used.
This will lead to the computation of 126 coefficients.

```
totalDegree = 5
enum_func = multivariateBasis.getEnumerateFunction()
P = enum_func.getBasisSizeFromTotalDegree(totalDegree)
print(f"P={P}")
```

```
P=126
```

We select the `FixedStrategy`

truncation rule, which corresponds to using
the first P polynomials of the polynomial basis.

```
adaptiveStrategy = ot.FixedStrategy(multivariateBasis, P)
```

We begin by getting the standard measure associated with the multivariate polynomial basis. We see that the range of the Beta distribution has been standardized into the [-1,1] interval. This is the same for the Uniform distribution and the second Beta distribution.

```
measure = multivariateBasis.getMeasure()
print(f"measure={measure}")
```

```
measure=JointDistribution(Beta(alpha = 0.9, beta = 3.5, a = -1, b = 1), LogNormal(muLog = 5.69881, sigmaLog = 0.0997513, gamma = 0), Uniform(a = -1, b = 1), Beta(alpha = 2.5, beta = 4, a = -1, b = 1), IndependentCopula(dimension = 4))
```

The choice of the `GaussProductExperiment`

rule with 4 nodes
in each of the 4 dimensions leads to evaluations of the model.

```
marginalSizes = [4] * dim_input
experiment = ot.GaussProductExperiment(distribution, marginalSizes)
print(f"N={experiment.getSize()}")
X, W = experiment.generateWithWeights()
Y = g(X)
```

```
N=256
```

We now set the method used to compute the coefficients; we select the integration method.

```
projectionStrategy = ot.IntegrationStrategy()
```

We can now create the functional chaos.

```
algo = ot.FunctionalChaosAlgorithm(
X, W, Y, distribution, adaptiveStrategy, projectionStrategy
)
algo.run()
```

Get the result

```
result = algo.getResult()
```

The `getMetaModel()`

method returns the metamodel function.

```
metamodel = result.getMetaModel()
```

## Validate the metamodel¶

Generate a new validation sample (which is independent of the training sample).

```
n_valid = 1000
X_test = distribution.getSample(n_valid)
Y_test = g(X_test)
```

The `MetaModelValidation`

class validates the metamodel
based on a validation sample.

```
val = ot.MetaModelValidation(Y_test, metamodel(X_test))
```

Compute the predictivity coefficient.

```
R2 = val.computeR2Score()[0]
R2
```

```
0.9999979360038378
```

Plot the observed versus the predicted outputs.

```
graph = val.drawValidation()
graph.setTitle(f"Gauss product N={experiment.getSize()} - R2={R2 * 100:.2f}")
view = otv.View(graph)
```

Now repeat the same process on various designs.

```
def draw_validation(experiment):
projectionStrategy = ot.IntegrationStrategy(experiment)
algo = ot.FunctionalChaosAlgorithm(
X, Y, distribution, adaptiveStrategy, projectionStrategy
)
algo.run()
result = algo.getResult()
metamodel = result.getMetaModel()
X_test = distribution.getSample(n_valid)
Y_test = g(X_test)
val = ot.MetaModelValidation(Y_test, metamodel(X_test))
R2 = val.computeR2Score()[0]
graph = val.drawValidation()
graph.setTitle(
f"{experiment.__class__.__name__} - N={experiment.getSize()} - R2={R2 * 100:.2f}"
)
return graph
```

Use an LHS design.

```
experiment = ot.LHSExperiment(distribution, int(1e6))
graph = draw_validation(experiment)
view = otv.View(graph)
```

Use a low-discrepancy experiment (Quasi-Monte Carlo).

```
sequence = ot.SobolSequence()
experiment = ot.LowDiscrepancyExperiment(sequence, distribution, int(1e5))
graph = draw_validation(experiment)
view = otv.View(graph)
otv.View.ShowAll()
```

## Conclusion¶

With the Gauss product rule the coefficients are particularly well computed since the Q2 coefficient is excellent, even with the relatively limited amount of simulation (256 points). On the other hand the LHS and low-discrepancy experiments require many more points to achieve a Q2>99%.