Calibration of the deflection of a tube

We consider a calibration of the deflection of a tube as described here.

from openturns.usecases import deflection_tube as deflection_tube
import openturns as ot
import openturns.viewer as viewer
from matplotlib import pylab as plt
ot.Log.Show(ot.Log.NONE)

Create a calibration problem

We load the model from the use case module :

dt = deflection_tube.DeflectionTube()

We create a sample out of our input distribution :

sampleSize = 100
inputSample = dt.inputDistribution.getSample(sampleSize)
inputSample[0:5]

	Force	Length	Location	External diameter	Internal diameter	Young Modulus
0	1.04272	1.5	1	0.8	0.1	200261.8
1	0.9817007	1.5	1	0.8	0.1	199645.5
2	1.030777	1.5	1	0.8	0.1	200796.7
3	1.001165	1.5	1	0.8	0.1	199523.4
4	0.9948547	1.5	1	0.8	0.1	197404.2

We take the image of our input sample by the model :

outputDeflection = dt.model(inputSample)
outputDeflection[0:5]

	Deflection	Left angle	Right angle
0	-7.195208e-06	-1.439042e-05	1.798802e-05
1	-6.795061e-06	-1.359012e-05	1.698765e-05
2	-7.093843e-06	-1.418769e-05	1.773461e-05
3	-6.934025e-06	-1.386805e-05	1.733506e-05
4	-6.964293e-06	-1.392859e-05	1.741073e-05

observationNoiseSigma = [0.1e-6, 0.05e-5, 0.05e-5]
observationNoiseCovariance = ot.CovarianceMatrix(3)
for i in range(3):
    observationNoiseCovariance[i, i] = observationNoiseSigma[i]**2

noiseSigma = ot.Normal([0., 0., 0.], observationNoiseCovariance)
sampleObservationNoise = noiseSigma.getSample(sampleSize)
observedOutput = outputDeflection + sampleObservationNoise
observedOutput[0:5]

	Deflection	Left angle	Right angle
0	-7.353985e-06	-1.414153e-05	1.7227e-05
1	-6.542875e-06	-1.434351e-05	1.789852e-05
2	-7.156474e-06	-1.452301e-05	1.771004e-05
3	-6.957683e-06	-1.433306e-05	1.735952e-05
4	-6.962969e-06	-1.336576e-05	1.802425e-05

observedInput = ot.Sample(sampleSize, 2)
observedInput[:, 0] = inputSample[:, 0]  # F
observedInput[:, 1] = inputSample[:, 5]  # E
observedInput.setDescription(["Force", "Young Modulus"])
observedInput[0:5]

	Force	Young Modulus
0	1.04272	200261.8
1	0.9817007	199645.5
2	1.030777	200796.7
3	1.001165	199523.4
4	0.9948547	197404.2

fullSample = ot.Sample(sampleSize, 5)
fullSample[:, 0:2] = observedInput
fullSample[:, 2:5] = observedOutput
fullSample.setDescription(
    ["Force", "Young", "Deflection", "Left Angle", "Right Angle"])
fullSample[0:5]

	Force	Young	Deflection	Left Angle	Right Angle
0	1.04272	200261.8	-7.353985e-06	-1.414153e-05	1.7227e-05
1	0.9817007	199645.5	-6.542875e-06	-1.434351e-05	1.789852e-05
2	1.030777	200796.7	-7.156474e-06	-1.452301e-05	1.771004e-05
3	1.001165	199523.4	-6.957683e-06	-1.433306e-05	1.735952e-05
4	0.9948547	197404.2	-6.962969e-06	-1.336576e-05	1.802425e-05

graph = ot.VisualTest.DrawPairs(fullSample)
view = viewer.View(graph)

Setting up the calibration

XL = 1.4  # Exact : 1.5
Xa = 1.2  # Exact : 1.0
XD = 0.7  # Exact : 0.8
Xd = 0.2  # Exact : 0.1
thetaPrior = [XL, Xa, XD, Xd]

sigmaXL = 0.1 * XL
sigmaXa = 0.1 * Xa
sigmaXD = 0.1 * XD
sigmaXd = 0.1 * Xd
parameterCovariance = ot.CovarianceMatrix(4)
parameterCovariance[0, 0] = sigmaXL**2
parameterCovariance[1, 1] = sigmaXa**2
parameterCovariance[2, 2] = sigmaXD**2
parameterCovariance[3, 3] = sigmaXd**2
parameterCovariance

[[ 0.0196 0 0 0 ]
[ 0 0.0144 0 0 ]
[ 0 0 0.0049 0 ]
[ 0 0 0 0.0004 ]]

calibratedIndices = [1, 2, 3, 4]
calibrationFunction = ot.ParametricFunction(
    dt.model, calibratedIndices, thetaPrior)

sigmaObservation = [0.2e-6, 0.03e-5, 0.03e-5]  # Exact : 0.1e-6

errorCovariance = ot.CovarianceMatrix(3)
errorCovariance[0, 0] = sigmaObservation[0]**2
errorCovariance[1, 1] = sigmaObservation[1]**2
errorCovariance[2, 2] = sigmaObservation[2]**2

calibrationFunction.setParameter(thetaPrior)
predictedOutput = calibrationFunction(observedInput)
predictedOutput[0:5]

	Deflection	Left angle	Right angle
0	-3.049682e-06	-1.016561e-05	1.651911e-05
1	-2.88008e-06	-9.600267e-06	1.560043e-05
2	-3.006719e-06	-1.00224e-05	1.628639e-05
3	-2.93898e-06	-9.7966e-06	1.591947e-05
4	-2.951809e-06	-9.839363e-06	1.598896e-05

Calibration with gaussian non linear least squares

algo = ot.GaussianNonLinearCalibration(
    calibrationFunction, observedInput, observedOutput, thetaPrior, parameterCovariance, errorCovariance)

algo.run()

calibrationResult = algo.getResult()

Analysis of the results

thetaMAP = calibrationResult.getParameterMAP()
thetaMAP

[1.49659,0.994631,0.800523,0.199881]

Compute a 95% confidence interval for each marginal.

thetaPosterior = calibrationResult.getParameterPosterior()
alpha = 0.95
dim = thetaPosterior.getDimension()
for i in range(dim):
    print(thetaPosterior.getMarginal(i).computeBilateralConfidenceInterval(alpha))

Out:

[1.47397, 1.51774]
[0.968328, 1.01897]
[0.795563, 0.805443]
[0.199878, 0.19991]

graph = calibrationResult.drawObservationsVsInputs()
view = viewer.View(graph)

graph = calibrationResult.drawObservationsVsPredictions()
view = viewer.View(graph)

graph = calibrationResult.drawResiduals()
view = viewer.View(graph)

graph = calibrationResult.drawParameterDistributions()
view = viewer.View(graph)

plt.show()