Calibration of the deflection of a tube

We calibrate the deflection of a tube as described here. More precisely, we calibrate the mechanical parameters of a model which computes the vertical deflection of a tube and two deflection angles. This example shows how to calibrate a computer code which has several outputs.

from openturns.usecases import 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	0.9886198	1.5	1	0.8	0.1	198813.2
1	0.9659264	1.5	1	0.8	0.1	197055.1
2	1.097242	1.5	1	0.8	0.1	199363.4
3	1.00796	1.5	1	0.8	0.1	199045.5
4	0.9531129	1.5	1	0.8	0.1	198060.4

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

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

	Deflection	Left angle	Right angle
0	-6.8716e-06	-1.37432e-05	1.7179e-05
1	-6.773763e-06	-1.354753e-05	1.693441e-05
2	-7.605552e-06	-1.52111e-05	1.901388e-05
3	-6.99785e-06	-1.39957e-05	1.749462e-05
4	-6.64998e-06	-1.329996e-05	1.662495e-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.0, 0.0], observationNoiseCovariance)
sampleObservationNoise = noiseSigma.getSample(sampleSize)
observedOutput = outputDeflection + sampleObservationNoise
observedOutput[0:5]

	Deflection	Left angle	Right angle
0	-6.770587e-06	-1.346099e-05	1.738157e-05
1	-6.710988e-06	-1.336128e-05	1.663024e-05
2	-7.714272e-06	-1.408485e-05	1.894873e-05
3	-6.932516e-06	-1.332838e-05	1.861578e-05
4	-6.626314e-06	-1.203596e-05	1.667933e-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	0.9886198	198813.2
1	0.9659264	197055.1
2	1.097242	199363.4
3	1.00796	199045.5
4	0.9531129	198060.4

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	0.9886198	198813.2	-6.770587e-06	-1.346099e-05	1.738157e-05
1	0.9659264	197055.1	-6.710988e-06	-1.336128e-05	1.663024e-05
2	1.097242	199363.4	-7.714272e-06	-1.408485e-05	1.894873e-05
3	1.00796	199045.5	-6.932516e-06	-1.332838e-05	1.861578e-05
4	0.9531129	198060.4	-6.626314e-06	-1.203596e-05	1.667933e-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	-2.912521e-06	-9.708404e-06	1.577616e-05
1	-2.871053e-06	-9.570177e-06	1.555154e-05
2	-3.223606e-06	-1.074535e-05	1.74612e-05
3	-2.966032e-06	-9.886773e-06	1.606601e-05
4	-2.818588e-06	-9.395293e-06	1.526735e-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.51028,1.01202,0.802726,0.199919]

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))

[1.49265, 1.52767]
[0.989565, 1.0335]
[0.799133, 0.80639]
[0.199886, 0.19997]

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()