# Calibration of an Unstable Dynamical System using Prediction Error Method

In this tutorial we show how can one use Prediction Error Method for calibration of an unstable dynamical system, namely the "ball and beam" system.

As shown in the schematic below, in this case, our process model involves, a horizontal beam and a motor. The control objective is to balance a rolling ball on the beam, when the position of the ball changes, by controlling the beam angle using a motor. Here, the system has three parameters namely ${I, g, f_v}$, where `I`

stands for the non-dimensional parameter, ${J/(m*r^2)}$, which is analytically 0.4 in this case for the spherical-ball, `g`

stands for the acceleration due to gravity, and ${f_v}$ stands for the Coulomb's friction coefficient parameter, with all units taken in S.I. Assuming there is negligible dynamics inside the motor components, hence the angular velocity of beam, is directly proportional to the input. The dynamical equations of the ball, based on the classical mechanics for a Coulomb's friction model, can be derived from the conservation equations described below.

The dataset for this system, has been collected from the ball-and-beam experiment setup in the STADIUS's Identification Database.

The measured signals out of this process, are two of the state variables of the problem, the beam angle relative to the horizontal plane (φ) and the position of the ball (x). First, we start the problem, by defining the governing system of ODEs, that consists of the following equations:-

\[\begin{aligned} \frac{d\phi}{dt} &= f^n(t) \\ (1+ I) \frac{d^2x}{dt^2} - x (\frac{d\phi}{dt})^2 + f_v \frac{dx}{dt} + g sin(\phi) &= 0 \\ \end{aligned}\]

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

ModelingToolkitStandardLibrary | Library for using standard modeling components |

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 |

DataInterpolations | Library for creating interpolations from the data |

ForwardDiff | Library for Forward Mode Automatic Differentiation |

Plots | The plotting and visualization library |

```
using JuliaSimModelOptimizer
using ModelingToolkit
import ModelingToolkit: D_nounits as D, t_nounits as t
using ModelingToolkitStandardLibrary.Blocks: RealInput, TimeVaryingFunction
using OrdinaryDiffEq
using CSV, DataFrames
using DataSets
using DataInterpolations: CubicSpline
using ForwardDiff: derivative
using Plots
```

## Data Setup

```
training_dataset = dataset("juliasimtutorials/ball_beam_data")
data = open(IO, training_dataset) do io
CSV.read(io,
DataFrame,
header = ["ballandbeamsystem.ϕ", "ballandbeamsystem.x"],
delim = ' ',
ignorerepeated = true)
end
df = hcat(DataFrame("timestamp" => range(0, 100, 1000)), data)
first(df, 5)
```

Row | timestamp | ballandbeamsystem.ϕ | ballandbeamsystem.x |
---|---|---|---|

Float64 | Float64 | Float64 | |

1 | 0.0 | -0.0019635 | -0.0048852 |

2 | 0.1001 | -0.0019635 | -0.00464094 |

3 | 0.2002 | -0.0019635 | -0.00378603 |

4 | 0.3003 | -0.0019635 | -0.00403029 |

5 | 0.4004 | -0.0019635 | -0.00390816 |

## Model Setup

```
t_vec = df[!, "timestamp"]
ϕ_vec = data[:, 1]
input_func = CubicSpline(ϕ_vec, t_vec)
dinput = t_vec -> derivative(input_func, t_vec)
dinput_vec = dinput.(t_vec)
dinput_func = CubicSpline(dinput_vec, t_vec)
pos_func = CubicSpline(data[:, 2], t_vec)
dpos_func = t_vec -> derivative(pos_func, t_vec)
dpos_vec = dpos_func.(t_vec)
function ballandbeamsystem()
states =@variables x(t)=-4.8851980e-03 ϕ(t)=-1.9634990e-03 ball_position(t)
ps = @parameters I=0.2 g=9.8 Fv=1.0
@named input = RealInput()
eqs = [
D(ϕ) ~ input.u,
0 ~ (1 + I) * D(D(x)) - x * (D(ϕ))^2 - g * sin(ϕ) + Fv * D(x),
]
@named ballandbeamsystem = ODESystem(eqs, t, states, ps, systems = [input])
end
function SystemModel()
f = dinput_func
bbsys = ballandbeamsystem()
t = ModelingToolkit.get_iv(bbsys)
src = TimeVaryingFunction(f; name = :src)
eqs = [ModelingToolkit.connect(bbsys.input, src.output)]
@named model = ODESystem(eqs, t; systems = [bbsys, src])
end
model = complete(SystemModel())
sys = structural_simplify(model)
```

The data that we have for the system corresponds to the angle of the beam relative to the horizontal plane (φ) and the position of the ball on the beam (x). Also, in this system, we need to pass the un-initialised dummy derivatives terms to the existing initialized variable map.

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

`experiment = Experiment(df, sys, alg = nothing, initial_conditions = [D(model.ballandbeamsystem.x) => 0])`

```
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, 100.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 = InverseProblem(experiment, [sys.ballandbeamsystem.I => (0.0, 10.0), sys.ballandbeamsystem.Fv => (0.0, 2.0)])`

```
InverseProblem with one experiment with 2 elements in the search space.
```

## Calibration

### SingleShooting

Let us first try to solve this problem using `SingleShooting`

. To do this, we first define an algorithm `alg`

and then call `calibrate`

with the `prob`

and `alg`

.

```
alg = SingleShooting(maxiters = 1000)
r = calibrate(prob, alg)
```

```
Calibration result computed in 2 minutes and 48 seconds. Final objective value: 4.14636.
Optimization ended with MaxIters Return Code and returned.
┌─────────────────────┬──────────────────────┐
│ ballandbeamsystem₊I │ ballandbeamsystem₊Fv │
├─────────────────────┼──────────────────────┤
│ 5.83433 │ 1.86311 │
└─────────────────────┴──────────────────────┘
```

The calibrated parameters don't look right. The value of `I`

and `Fv`

seems to be very high. We can use these parameters to simulate it and plot to see if it fits the data well.

`plot(experiment, prob, r, show_data = true, ms = 1.0, size = (1000, 600), layout = (2, 1))`

### SingleShooting with Prediction Error Method

We can see that the simulation does not match the data at all! It is because this is an unstable system and hence simulating is very difficult. So, the calibration process is also very difficult as simulations will diverge and we can never fit the data correctly. To mitigate this, we will use Prediction Error Method where the simulation is guided by the data such that the trajectory won't diverge and this should help with the calibration process.

To use Prediction Error Method, we need to pass it in the `model_transformations`

keyword in the constructor of the `Experiment`

.

```
experiment = Experiment(df, sys, alg = nothing, initial_conditions = [D(model.ballandbeamsystem.x) => 0], model_transformations = [DiscreteFixedGainPEM(0.2)])
prob = InverseProblem(experiment, [sys.ballandbeamsystem.I => (0.0, 10.0), sys.ballandbeamsystem.Fv => (0.0, 2.0)])
```

```
InverseProblem with one experiment with 2 elements in the search space.
```

Argument passed to `DiscreteFixedGainPEM`

is the amount of correction needed during simulation. `1.0`

represents completely using the data and `0.0`

represents completely ignoring the data. Typically, we should use this be about 0.2-0.3 to help guide the simulation.

Now, if we try calibrating, we get:

```
alg = SingleShooting(maxiters = 1000)
r = calibrate(prob, alg)
```

```
Calibration result computed in 44 seconds. Final objective value: 0.037686.
Optimization ended with Success Return Code and returned.
┌─────────────────────┬──────────────────────┐
│ ballandbeamsystem₊I │ ballandbeamsystem₊Fv │
├─────────────────────┼──────────────────────┤
│ 2.07774 │ 0.722475 │
└─────────────────────┴──────────────────────┘
```

The calibrated value of `I`

looks a bit high but it is closer to its theoretical value from using only `SingleShooting`

. We can again use the calibrated parameters to simulate it and see if it fits the data well.

`plot(experiment, prob, r, show_data = true, ms = 1.0, size = (1000, 600), layout = (2, 1))`

It does fit the data well! The calibrated parameters look a bit high than their theoretical values as the model may not be a sufficient representation of the real world experiment from where the data was collected. One way to improve the results to have a more sophisticated model to improve the results, i.e., get the values of the calibrated parameters closer to what they actually represent. But one can see that using Prediction Error Method, the results improve drastically, from not fitting the data at all to closely fitting the data which also helps in the calibration process.