Calibrate a parametric model: a quick-start guide to calibration

In this example we present the calibration of a parametric module. To do this, we show how to define the observations. Then we create a Python function and create the parametric model that is to be calibrated. Finally, we calibrate the parameters of the model using least squares and we validate the hypotheses of the method.

Please read Code calibration for more details on code calibration and least squares. In this example we are interested in the calibration of the flooding model. Once the reader has mastered this example, the calibration of the Chaboche mechanical model may be considered to make an in-depth study of these algorithms. The goal of this script is to be relatively small, so please consider reading the other examples if this is relevant.

Parameters to calibrate and observations

The parametric model is:

The vector of parameters to calibrate is:

We consider a sample size equal to:

The observations are the couples , i.e. each observation is a couple made of the flowrate and the corresponding river height.

We choose to calibrate this model using non linear least squares, because this is a relatively flexible method.

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

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

Define the observations

In practice, we generally use a data set which has been obtained from measurements. This data set can be loaded using e.g. ImportFromCSVFile(). Here we import the data from the FloodModel class.

fm = flood_model.FloodModel()
Qobs = fm.data[:, 0]
Hobs = fm.data[:, 1]
nbobs = fm.data.getSize()
fm.data

	Q ($m^3/s$)	H (m)
0	130	0.59
1	530	1.33
2	960	2.03
3	1400	2.72
4	1830	2.83
5	2260	3.5
6	2700	3.82
7	3130	4.36
8	3560	4.63
9	4010	4.96

Create the physical model

We define the model which has 3 inputs and one output H. This model has a parametric shape that may correspond to the data, which has some power shape. In the model, the parameters are a and b, the input is Q and the output is H: - Q : the flowrate of the river, - a, b : the parameters.

def functionSimpleFlooding(X):
    Q, a, b = X
    H = a * Q**b
    return [H]


g = ot.PythonFunction(3, 1, functionSimpleFlooding)
g = ot.MemoizeFunction(g)
g.setInputDescription(["Q", "a", "b"])
g.setOutputDescription(["H"])

Setting the calibration parameters

Define the value of the reference values of the parameter. There is no particular method to set these values: we used trial-and-error to see the order of magnitude of the parameters.

aInitial = 0.1
bInitial = 0.5
thetaPrior = [aInitial, bInitial]

Create the parametric function

In the physical model, the inputs and parameters are ordered as presented in the next table. Notice that there are no parameters in the physical model.

Index	Input variable
0	Q
1	a
2	b

Index	Parameter
∅	∅

Table 1. Indices and names of the inputs and parameters of the physical model.

print("Physical Model Inputs:", g.getInputDescription())
print("Physical Model Parameters:", g.getParameterDescription())

Physical Model Inputs: [Q,a,b]
Physical Model Parameters: []

In order to perform calibration, we have to define a parametric model, with observed inputs and parameters to calibrate. In order to do this, we create a ParametricFunction where the parameters are a and b which have the indices 1 and and 2 in the physical model. Please read Create a parametric function for more details on this topic.

Index	Input variable
0	Q

Index	Parameter
0	a
1	b

Table 2. Indices and names of the inputs and parameters of the parametric model.

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

calibratedIndices = [1, 2]
mycf = ot.ParametricFunction(g, calibratedIndices, thetaPrior)

Plot the Y observations versus the X observations.

title = "Before calibration : a = %.4f, b = %.4f" % (aInitial, bInitial)
graph = ot.Graph(title, "Q", "H", True)
# Plot the model before calibration
curve = mycf.draw(100.0, 4000.0).getDrawable(0)
curve.setLegend("Model, before calibration")
graph.add(curve)
# Plot the noisy observations
cloud = ot.Cloud(Qobs, Hobs)
cloud.setLegend("Observations")
graph.add(cloud)
#
graph.setLegendPosition("upper left")
view = otv.View(graph)

Wee see that the model does not fit to the data. The goal of calibration is to find which parameter best fit to the observations.

Calibration with non linear least squares

The NonLinearLeastSquaresCalibration class performs non linear least squares calibration by minimizing the squared Euclidian norm between the predictions and the observations.

algo = ot.NonLinearLeastSquaresCalibration(mycf, Qobs, Hobs, thetaPrior)

The run() method computes the solution of the problem.

algo.run()
calibrationResult = algo.getResult()

Analysis of the results

The getParameterMAP() method returns the maximum of the posterior distribution of .

thetaMAP = calibrationResult.getParameterMAP()
print(thetaMAP)

[0.0258034,0.634465]

In order to see if the parameters fit the data, we plot the input depending on the output before and after calibration.

# sphinx_gallery_thumbnail_number = 2
graph = calibrationResult.drawObservationsVsInputs()
aEstimated, bEstimated = thetaMAP
title = "After calibration : a = %.4f, b = %.4f" % (aEstimated, bEstimated)
graph.setTitle(title)
graph.setLegendPosition("upper left")
view = otv.View(graph)

One of the hypotheses of the least squares method is that the residuals follow a normal distribution: the next cell checks if this hypothesis is satisfied here.

graph = calibrationResult.drawResiduals()
graph.setLegendPosition("upper left")
view = otv.View(graph)

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.

Remark. The logarithm of the height is:

Hence, transforming the data on a logarithmic scale leads to a parametric model that is linear in the parameters. The parameters of this transformed model can be estimated using linear linear squares, which may lead to a significant improvement in terms of number of function evaluations.

otv.View.ShowAll()