# 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:

Module | Description |
---|---|

JuliaSimModelOptimizer | The high-level library used to formulate our problem and perform automated model discovery |

ModelingToolkit | The symbolic modeling environment |

OrdinaryDiffEq | The numerical differential equation solvers |

CSV and DataFrames | We will read our experimental data from .csv files |

DataSets | We will load our experimental data from datasets on JuliaHub |

Plots | The plotting and visualization library |

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

Row | timestamp | m |
---|---|---|

Float64 | Float64 | |

1 | 0.1 | 0.828499 |

2 | 0.11992 | 0.957478 |

3 | 0.13984 | 1.06817 |

4 | 0.15976 | 1.15853 |

5 | 0.17968 | 1.22694 |

## 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 48 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 810 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.