Calibration of a Dynamical System with Observed Data with Prediction Error Method

In this tutorial, we will show how to calibrate length of a pendulum using prediction error method.

Prediction error methods (PEM) are good at smoothing loss landscapes which eases the calibration process. Without PEM, this problem has tons of local minima in the loss landscape which makes the calibration process difficult to converge to the correct value.

Julia Environment

For this tutorial we will need the following packages:

using JuliaSimModelOptimizer
using ModelingToolkit
import ModelingToolkit: D_nounits as D, t_nounits as t
using OrdinaryDiffEq
using CSV, DataFrames
using DataSets
using Plots

Data Setup

Let us import the dataset. This dataset contains data for the observed variables of the system.

pendulum_dataset = dataset("juliasimtutorials/pendulum_data")

data = open(IO, pendulum_dataset) do io
    CSV.read(io, DataFrame)
end

first(data, 5)

Model Setup

We will now define the model using ModelingToolkit. The model is a simple pendulum fixed at one end and mass on the other. We will also assume mass of the string is negiglible. x₁ denotes the displacement and x₂ denotes the velocity. We are measuring 3*x₁ which is the data collected.

@variables x₁(t)=0.0 x₂(t)=3.0 m(t)
@parameters L = 2.0
@constants g = 9.81
tspan = (0.0, 20.0)
eqs = [
    D(x₁) ~ x₂,
    D(x₂) ~ -(g / L) * sin(x₁),
    m ~ 3 * x₁
]
@named model = ODESystem(eqs, t, [x₁, x₂], [L]; tspan)
sys = structural_simplify(model)

\[ \begin{align} \frac{\mathrm{d} x_1\left( t \right)}{\mathrm{d}t} =& x_2\left( t \right) \\ \frac{\mathrm{d} x_2\left( t \right)}{\mathrm{d}t} =& \frac{ - g \sin\left( x_1\left( t \right) \right)}{L} \end{align} \]

Defining Experiment and InverseProblem

In order to create an Experiment, we will use the default initial values of the states and parameters of our model. These are our initial guesses which will be used to optimize the inverse problem in order to fit the given data. To use prediction error method, we pass DiscreteFixedGainPEM in the model_transformations keyword argument. DiscreteFixedGainPEM automatically tries to use all the data present for guiding the simulation - both unknowns and observed. If observed data is present, it tries to solve the system of equations to estimate unknowns, provided the observed equations corresponding to the data are linear and invertible.

experiment_no_pem = Experiment(data, sys; abstol = 1e-8, reltol = 1e-6)
experiment_pem = Experiment(data, sys; reltol = 1e-6, model_transformations = [DiscreteFixedGainPEM(0.1)])

Experiment for model with no fixed parameters or initial conditions.
The simulation of this experiment is given by:
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 20.0)

Once we have created the experiment, the next step is to create an InverseProblem. This inverse problem, requires us to provide the search space as a vector of pairs corresponding to the parameters that we want to recover and the assumption that we have for their respective bounds.

prob_no_pem = InverseProblem(experiment_no_pem, [L => (0.0, 3.0)])
prob_pem = InverseProblem(experiment_pem, [L => (0.0, 3.0)])

InverseProblem with one experiment with 1 elements in the search space.

Calibration

The true length of the pendulum is 0.2.

SingleShooting

Let us first try to solve this problem using SingleShooting without PEM. To do this, we first define an algorithm alg and then call calibrate with the prob_no_pem and alg.

alg = SingleShooting(maxiters = 1000)
r_no_pem = calibrate(prob_no_pem, alg)

Calibration result computed in 52 seconds. Final objective value: 841.545.
Optimization ended with Success Return Code and returned.

┌────────────┐
│          L │
├────────────┤
│ 7.86704e-7 │
└────────────┘

We can see the calibrated length of the pendulum is far from the correct value.

SingleShooting with Prediction Error Method

Now, let us try it with PEM.

alg = SingleShooting(maxiters = 1000)
r_pem = calibrate(prob_pem, alg)

Calibration result computed in 29 seconds and 284 milliseconds. Final objective value: 4.53701e-12.
Optimization ended with Success Return Code and returned.

┌─────┐
│   L │
├─────┤
│ 0.2 │
└─────┘

We can see it calibrates correctly!

Visualization

Let us plot the result to confirm visually.

plot(experiment_no_pem, prob_no_pem, r_no_pem, show_data = true, ms = 1.0, size = (1000, 300))

plot(experiment_pem, prob_pem, r_pem, show_data = true, ms = 1.0, size = (1000, 300))

We can see the simulation results without PEM do not match the data at all but with PEM matches the data really well.