Generate observations of the Chaboche mechanical model

In this example we present the simulation of noisy observations of the Chaboche model. A detailed explanation of this mechanical law is presented here. We show how to produce the observations that we use in the calibration model of Calibration of the Chaboche model.

Parameters to calibrate

The vector of parameters to calibrate is:

\theta = (R,C,\gamma).

The true values of the parameters are:

  • R = 750\times 10^6,

  • C = 2750\times 10^6,

  • \gamma = 10.

This is the set of true values that we wish to estimate with the calibration methods. In practical studies, these values are unknown.

Observations

In order to create a calibration problem, we make the hypothesis that the strain has the following distribution:

\epsilon \sim \mathcal{U}(0,0.07).

where \mathcal{U} is the uniform distribution. Moreover, we consider a Gaussian noise on the observed constraint:

\epsilon_\sigma \sim \mathcal{N} \left(0,10\times 10^6\right)

and we make the hypothesis that the observation errors are independent. We set the number of observations to:

n = 10.

We generate a Monte-Carlo sample with size n:

\sigma_i = g(\epsilon_i,R,C,\gamma) + (\epsilon_\sigma)_i,

for i = 1,..., n. The observations are the pairs \{(\epsilon_i,\sigma_i)\}_{i=1,...,n}, i.e. each observation is a couple made of the strain and the corresponding stress.

Variables

In the particular situation where we want to calibrate this model, the following list presents which variables are observed input variables, input calibrated variables and observed output variables.

  • \epsilon: Input. Observed.

  • R, C, \gamma : Inputs. Calibrated.

  • \sigma: Output. Observed.

import openturns as ot
import openturns.viewer as otv
from openturns.usecases import chaboche_model

ot.Log.Show(ot.Log.NONE)

Generate the observations

In practice, we generally use a data set which has been obtained from measurements. In this example, we generate the data using noisy observations of the physical model. In the next part, we will calibrate the parameters using the calibration algorithms.

Load the Chaboche model

cm = chaboche_model.ChabocheModel()
print("Inputs:", cm.model.getInputDescription())
print("Outputs:", cm.model.getOutputDescription())
Inputs: [Strain,R,C,Gamma]
Outputs: [Stress]

We see that there are four inputs: Strain, R, C and Gamma and one output: Sigma. The Strain is observed on input and the stress Sigma is observed on output. Using these observations, we want to estimate the parameters R, C and Gamma.

We get the Chaboche model and the joint input distribution :

inputDistribution = cm.inputDistribution
print("inputDistribution:")
inputDistribution
inputDistribution:
JointDistribution
  • name=JointDistribution
  • dimension: 4
  • description=[Strain,R,C,Gamma]
  • copula: IndependentCopula(dimension = 4)
Index Variable Distribution
0 Strain Uniform(a = 0, b = 0.07)
1 R LogNormal(muLog = 20.4355, sigmaLog = 0.0146659, gamma = 0)
2 C Normal(mu = 2.75e+09, sigma = 2.5e+08)
3 Gamma Normal(mu = 10, sigma = 2)


Set the calibration parameters

In this part, we begin the calibration study. Define the value of the reference values of the \theta parameter. In the Bayesian framework, this is called the mean of the prior Gaussian distribution. In the data assimilation framework, this is called the background.

thetaTrue = [cm.trueR, cm.trueC, cm.trueGamma]
print("theta True = ")
print("    R = %.2f (MPa)" % (cm.trueR / 1.0e6))
print("    C = %.2f (MPa)" % (cm.trueC / 1.0e6))
print("    Gamma = %.4f" % (cm.trueGamma))
theta True =
    R = 750.00 (MPa)
    C = 2750.00 (MPa)
    Gamma = 10.0000

The following statement create the calibrated function from the model. The calibrated parameters R, C, Gamma are at indices 1, 2, 3 in the inputs arguments of the model.

calibratedIndices = [1, 2, 3]
mycf = ot.ParametricFunction(cm.model, calibratedIndices, thetaTrue)

Create a regular grid of the strains and evaluate the corresponding outputs.

nbobs = 10
step = (cm.strainMax - cm.strainMin) / (nbobs - 1)
rg = ot.RegularGrid(cm.strainMin, step, nbobs)
observedStrain = rg.getVertices()
predictedStress = mycf(observedStrain)

Generate observation noise. Here we consider a Gaussian observation noise, that we add to the output of the model.

stressObservationNoiseSigma = 10.0e6  # (Pa)
noiseSigma = ot.Normal(0.0, stressObservationNoiseSigma)
sampleNoiseStress = noiseSigma.getSample(nbobs)
observedStress = predictedStress + sampleNoiseStress

Gather the data into a sample.

data = ot.Sample(nbobs, 2)
data[:, 0] = observedStrain
data[:, 1] = observedStress
data.setDescription(["Strain", "Stress (Pa)"])
data
StrainStress (Pa)
007.568656e+08
10.0077777787.802e+08
20.015555567.945049e+08
30.023333337.976711e+08
40.031111118.357973e+08
50.038888898.181291e+08
60.046666678.40723e+08
70.054444448.668916e+08
80.062222228.68183e+08
90.078.927765e+08


Then we plot the model and compare it to the observations.

graph = ot.Graph("Model with true parameters", "Strain", "Stress (MPa)", True)
# Plot the model
curve = mycf.draw(cm.strainMin, cm.strainMax, 50).getDrawable(0)
curve.setLegend("True model")
curve.setLineStyle(ot.ResourceMap.GetAsString("CalibrationResult-ObservationLineStyle"))
graph.add(curve)
# Plot the noisy observations
cloud = ot.Cloud(observedStrain, observedStress)
cloud.setLegend("Observations")
cloud.setPointStyle(
    ot.ResourceMap.GetAsString("CalibrationResult-ObservationPointStyle")
)
graph.add(cloud)
graph.setLegendPosition("upper left")
view = otv.View(graph)
Model with true parameters

We see that the observations are relatively noisy, but that the trend is clear: this shows that it may be possible to fit the model. At this point, we have a data set that we can use for calibration and a model to calibrate.

otv.View.ShowAll()