`.
.. code-block:: default
import numpy as np
import openturns as ot
import openturns.viewer as viewer
from matplotlib import pylab as plt
from openturns.usecases import chaboche_model as chaboche_model
ot.Log.Show(ot.Log.NONE)
Load the Chaboche data structure
.. code-block:: default
cm = chaboche_model.ChabocheModel()
We get the Chaboche model and the joint input distribution :
.. code-block:: default
g = cm.model
inputDistribution = cm.inputDistribution
Create the Monte-Carlo sample.
.. code-block:: default
sampleSize = 100
inputSample = inputDistribution.getSample(sampleSize)
outputStress = g(inputSample)
outputStress[0:5]
.. raw:: html
| y0 |
|---|
| 0 | 8.225139e+08 |
|---|
| 1 | 7.87381e+08 |
|---|
| 2 | 7.530167e+08 |
|---|
| 3 | 8.25257e+08 |
|---|
| 4 | 8.215074e+08 |
|---|
Plot the histogram of the output.
.. code-block:: default
histoGraph = ot.HistogramFactory().build(outputStress / 1.0e6).drawPDF()
histoGraph.setTitle("Histogram of the sample stress")
histoGraph.setXTitle("Stress (MPa)")
histoGraph.setLegends([""])
view = viewer.View(histoGraph)
.. image:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_chaboche_001.png
:alt: Histogram of the sample stress
:class: sphx-glr-single-img
Generate observation noise.
.. code-block:: default
stressObservationNoiseSigma = 10.0e6 # (Pa)
noiseSigma = ot.Normal(0.0, stressObservationNoiseSigma)
sampleNoiseH = noiseSigma.getSample(sampleSize)
observedStress = outputStress + sampleNoiseH
.. code-block:: default
observedStrain = inputSample[:, 0]
.. code-block:: default
graph = ot.Graph("Observations", "Strain", "Stress (MPa)", True)
cloud = ot.Cloud(observedStrain, observedStress / 1.0e6)
graph.add(cloud)
view = viewer.View(graph)
.. image:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_chaboche_002.png
:alt: Observations
:class: sphx-glr-single-img
Set the calibration parameters
------------------------------
Define the value of the reference values of the :math:`\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*.
.. code-block:: default
R = 700e6 # Exact : 750e6
C = 2500e6 # Exact : 2750e6
Gamma = 8.0 # Exact : 10
thetaPrior = ot.Point([R, C, Gamma])
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.
.. code-block:: default
calibratedIndices = [1, 2, 3]
mycf = ot.ParametricFunction(g, calibratedIndices, thetaPrior)
Calibration with linear least squares
-------------------------------------
The `LinearLeastSquaresCalibration` class performs the linear least squares calibration by linearizing the model in the neighbourhood of the reference point.
.. code-block:: default
algo = ot.LinearLeastSquaresCalibration(
mycf, observedStrain, observedStress, thetaPrior, "SVD"
)
The `run` method computes the solution of the problem.
.. code-block:: default
algo.run()
.. code-block:: default
calibrationResult = algo.getResult()
Analysis of the results
-----------------------
The `getParameterMAP` method returns the maximum of the posterior distribution of :math:`\theta`.
.. code-block:: default
thetaMAP = calibrationResult.getParameterMAP()
thetaMAP
.. raw:: html
[7.55665e+08,2.41976e+09,7.76849]
We can compute a 95% confidence interval of the parameter :math:`\theta^\star`.
.. code-block:: default
thetaPosterior = calibrationResult.getParameterPosterior()
thetaPosterior.computeBilateralConfidenceIntervalWithMarginalProbability(0.95)[0]
.. raw:: html
[7.5e+08, 7.61329e+08]
[1.93278e+09, 2.90674e+09]
[-479.264, 494.801]
We can see that the `Gamma` parameter has a large confidence interval : even the sign of the parameter is unknown.
The parameter which is calibrated with the smallest confidence
interval is `R`, which confidence interval is [708.3,780.0] (MPa).
This is why this parameter seems the most important in this case.
.. code-block:: default
graph = calibrationResult.drawObservationsVsInputs()
graph.setLegendPosition("topleft")
view = viewer.View(graph)
.. image:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_chaboche_003.png
:alt: plot calibration chaboche
:class: sphx-glr-single-img
We see that there is a good fit after calibration, since the predictions after calibration (i.e. the green crosses) are close to the observations (i.e. the blue crosses).
.. code-block:: default
graph = calibrationResult.drawObservationsVsPredictions()
view = viewer.View(graph)
.. image:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_chaboche_004.png
:alt: plot calibration chaboche
:class: sphx-glr-single-img
We see that there is a much better fit after calibration, since the predictions are close to the diagonal of the graphics.
.. code-block:: default
observationError = calibrationResult.getObservationsError()
observationError
.. raw:: html
Normal(mu = 0, sigma = 1.03405e+07)
.. code-block:: default
graph = calibrationResult.drawResiduals()
graph.setLegendPosition("topleft")
view = viewer.View(graph)
.. image:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_chaboche_005.png
:alt: plot calibration chaboche
:class: sphx-glr-single-img
The analysis of the residuals shows that the distribution is centered on zero and symmetric. This indicates that the calibration performed well.
Moreover, the distribution of the residuals is close to being gaussian.
.. code-block:: default
graph = calibrationResult.drawParameterDistributions()
view = viewer.View(graph)
.. image:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_chaboche_006.png
:alt: plot calibration chaboche
:class: sphx-glr-single-img
Calibration with nonlinear least squares
----------------------------------------
The `NonLinearLeastSquaresCalibration` class performs the non linear least squares calibration by minimizing the squared euclidian norm between the predictions and the observations.
.. code-block:: default
algo = ot.NonLinearLeastSquaresCalibration(
mycf, observedStrain, observedStress, thetaPrior
)
The `run` method computes the solution of the problem.
.. code-block:: default
algo.run()
.. code-block:: default
calibrationResult = algo.getResult()
Analysis of the results
-----------------------
The `getParameterMAP` method returns the maximum of the posterior distribution of :math:`\theta`.
.. code-block:: default
thetaMAP = calibrationResult.getParameterMAP()
thetaMAP
.. raw:: html
[7.55764e+08,2.41907e+09,7.77538]
We can compute a 95% confidence interval of the parameter :math:`\theta^\star`.
.. code-block:: default
thetaPosterior = calibrationResult.getParameterPosterior()
thetaPosterior.computeBilateralConfidenceIntervalWithMarginalProbability(0.95)[0]
.. raw:: html
[7.46815e+08, 7.64023e+08]
[1.83849e+09, 3.1899e+09]
[0.976967, 16.4257]
We can see that all three parameters are estimated with a large confidence interval.
.. code-block:: default
graph = calibrationResult.drawObservationsVsInputs()
graph.setLegendPosition("topleft")
view = viewer.View(graph)
.. image:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_chaboche_007.png
:alt: plot calibration chaboche
:class: sphx-glr-single-img
We see that there is a good fit after calibration, since the predictions after calibration (i.e. the green crosses) are close to the observations (i.e. the blue crosses).
.. code-block:: default
graph = calibrationResult.drawObservationsVsPredictions()
view = viewer.View(graph)
.. image:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_chaboche_008.png
:alt: plot calibration chaboche
:class: sphx-glr-single-img
We see that there is a much better fit after calibration, since the predictions are close to the diagonal of the graphics.
.. code-block:: default
observationError = calibrationResult.getObservationsError()
observationError
.. raw:: html
Normal(mu = -18951.8, sigma = 1.0238e+07)
.. code-block:: default
graph = observationError.drawPDF()
view = viewer.View(graph)
.. image:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_chaboche_009.png
:alt: plot calibration chaboche
:class: sphx-glr-single-img
.. code-block:: default
graph = calibrationResult.drawResiduals()
graph.setLegendPosition("topleft")
view = viewer.View(graph)
.. image:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_chaboche_010.png
:alt: plot calibration chaboche
:class: sphx-glr-single-img
The analysis of the residuals shows that the distribution is centered on zero and symmetric. This indicates that the calibration performed well.
Moreover, the distribution of the residuals is close to being gaussian. This indicates that the observation error might be gaussian.
.. code-block:: default
graph = calibrationResult.drawParameterDistributions()
view = viewer.View(graph)
.. image:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_chaboche_011.png
:alt: plot calibration chaboche
:class: sphx-glr-single-img
Gaussian calibration parameters
-------------------------------
The standard deviation of the observations errors.
.. code-block:: default
sigmaStress = 1.0e7 # (Pa)
Define the covariance matrix of the output Y of the model.
.. code-block:: default
errorCovariance = ot.CovarianceMatrix(1)
errorCovariance[0, 0] = sigmaStress ** 2
Define the covariance matrix of the parameters :math:`\theta` to calibrate.
.. code-block:: default
sigmaR = 0.1 * R
sigmaC = 0.1 * C
sigmaGamma = 0.1 * Gamma
.. code-block:: default
sigma = ot.CovarianceMatrix(3)
sigma[0, 0] = sigmaR ** 2
sigma[1, 1] = sigmaC ** 2
sigma[2, 2] = sigmaGamma ** 2
sigma
.. raw:: html
[[ 4.9e+15 0 0 ]
[ 0 6.25e+16 0 ]
[ 0 0 0.64 ]]
Gaussian linear calibration
---------------------------
The `GaussianLinearCalibration` class performs the gaussian linear calibration by linearizing the model in the neighbourhood of the prior.
.. code-block:: default
algo = ot.GaussianLinearCalibration(
mycf, observedStrain, observedStress, thetaPrior, sigma, errorCovariance
)
The `run` method computes the solution of the problem.
.. code-block:: default
algo.run()
.. code-block:: default
calibrationResult = algo.getResult()
Analysis of the results
-----------------------
The `getParameterMAP` method returns the maximum of the posterior distribution of :math:`\theta`.
.. code-block:: default
thetaMAP = calibrationResult.getParameterMAP()
thetaMAP
.. raw:: html
[7.55403e+08,2.44237e+09,8.02337]
We can compute a 95% confidence interval of the parameter :math:`\theta^\star`.
.. code-block:: default
thetaPosterior = calibrationResult.getParameterPosterior()
thetaPosterior.computeBilateralConfidenceIntervalWithMarginalProbability(0.95)[0]
.. raw:: html
[7.51146e+08, 7.59661e+08]
[2.26218e+09, 2.62255e+09]
[6.28201, 9.76473]
We can see that all three parameters are estimated with a large confidence interval.
.. code-block:: default
graph = calibrationResult.drawObservationsVsInputs()
graph.setLegendPosition("topleft")
view = viewer.View(graph)
.. image:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_chaboche_012.png
:alt: plot calibration chaboche
:class: sphx-glr-single-img
We see that there is a good fit after calibration, since the predictions after calibration (i.e. the green crosses) are close to the observations (i.e. the blue crosses).
.. code-block:: default
graph = calibrationResult.drawObservationsVsPredictions()
view = viewer.View(graph)
.. image:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_chaboche_013.png
:alt: plot calibration chaboche
:class: sphx-glr-single-img
We see that there is a much better fit after calibration, since the predictions are close to the diagonal of the graphics.
The observation error is predicted by linearizing the problem at the prior.
.. code-block:: default
observationError = calibrationResult.getObservationsError()
observationError
.. raw:: html
Normal(mu = 0, sigma = 1e+07)
This can be compared to the residuals distribution, which is computed at the posterior.
.. code-block:: default
graph = calibrationResult.drawResiduals()
graph.setLegendPosition("topleft")
view = viewer.View(graph)
.. image:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_chaboche_014.png
:alt: plot calibration chaboche
:class: sphx-glr-single-img
The analysis of the gaussian distribution (the blue line) of the observation errors is close to the posterior distribution of the residuals (the green line). Moreover, the posterior distribution is centered. These informations indicate that the calibration performed well.
.. code-block:: default
graph = calibrationResult.drawParameterDistributions()
view = viewer.View(graph)
.. image:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_chaboche_015.png
:alt: plot calibration chaboche
:class: sphx-glr-single-img
For the :math:`R` variable, the observations are much more important than the prior: this is shown by the fact that the posterior and prior distribution of the :math:`R` variable are very different.
We see that the prior and posterior distribution are close to each other for the :math:`\gamma` variable: the observations did not convey much information for this variable.
Gaussian nonlinear calibration
------------------------------
The `GaussianNonLinearCalibration` class performs the gaussian nonlinear calibration.
.. code-block:: default
algo = ot.GaussianNonLinearCalibration(
mycf, observedStrain, observedStress, thetaPrior, sigma, errorCovariance
)
The `run` method computes the solution of the problem.
.. code-block:: default
algo.run()
.. code-block:: default
calibrationResult = algo.getResult()
Analysis of the results
-----------------------
The `getParameterMAP` method returns the maximum of the posterior distribution of :math:`\theta`.
.. code-block:: default
thetaMAP = calibrationResult.getParameterMAP()
thetaMAP
.. raw:: html
[7.5559e+08,2.43721e+09,8.00266]
We can compute a 95% confidence interval of the parameter :math:`\theta^\star`.
.. code-block:: default
thetaPosterior = calibrationResult.getParameterPosterior()
thetaPosterior.computeBilateralConfidenceIntervalWithMarginalProbability(0.95)[0]
.. raw:: html
[7.50579e+08, 7.60352e+08]
[2.2832e+09, 2.59508e+09]
[7.44971, 8.47513]
We can see that all three parameters are estimated with a large confidence interval.
.. code-block:: default
graph = calibrationResult.drawObservationsVsInputs()
graph.setLegendPosition("topleft")
view = viewer.View(graph)
.. image:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_chaboche_016.png
:alt: plot calibration chaboche
:class: sphx-glr-single-img
We see that there is a good fit after calibration, since the predictions after calibration (i.e. the green crosses) are close to the observations (i.e. the blue crosses).
.. code-block:: default
graph = calibrationResult.drawObservationsVsPredictions()
view = viewer.View(graph)
.. image:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_chaboche_017.png
:alt: plot calibration chaboche
:class: sphx-glr-single-img
We see that there is a much better fit after calibration, since the predictions are close to the diagonal of the graphics.
The observation error is predicted by bootstraping the problem at the posterior.
.. code-block:: default
observationError = calibrationResult.getObservationsError()
observationError
.. raw:: html
Normal(mu = -10638.6, sigma = 1e+07)
This can be compared to the residuals distribution, which is computed at the posterior.
.. code-block:: default
graph = calibrationResult.drawResiduals()
view = viewer.View(graph)
.. image:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_chaboche_018.png
:alt: plot calibration chaboche
:class: sphx-glr-single-img
The analysis is very similar to the linear calibration.
.. code-block:: default
graph = calibrationResult.drawParameterDistributions()
view = viewer.View(graph)
plt.show()
.. image:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_chaboche_019.png
:alt: plot calibration chaboche
:class: sphx-glr-single-img
We see that the prior and posterior distribution for the :math:`\gamma` parameter are close to each other, but not superimposed: the observations significantly brought information to the variable :math:`\gamma` during the calibration.
.. rst-class:: sphx-glr-timing
**Total running time of the script:** ( 0 minutes 8.521 seconds)
.. _sphx_glr_download_auto_calibration_least_squares_and_gaussian_calibration_plot_calibration_chaboche.py:
.. only :: html
.. container:: sphx-glr-footer
:class: sphx-glr-footer-example
.. container:: sphx-glr-download sphx-glr-download-python
:download:`Download Python source code: plot_calibration_chaboche.py `
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download Jupyter notebook: plot_calibration_chaboche.ipynb `
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery `_