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:

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)

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 1 minute and 27 seconds. Final objective value: 4.08715.
Optimization ended with MaxIters Return Code and returned.

┌─────────────────────┬──────────────────────┐
│ ballandbeamsystem₊I │ ballandbeamsystem₊Fv │
├─────────────────────┼──────────────────────┤
│             5.44539 │                  2.0 │
└─────────────────────┴──────────────────────┘

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