Parameter Estimation with Multiple Experiments

In this example, we will perform parameter calibration with two different experiments. The system we will calibrate is a double-mass system, where two masses are connected by a flexible shaft. An input force excites the first mass, which in turn drives the second mass through the flexible shaft. The flexibility of the shaft, modeled as a spring and damper, leads to a resonance, the properties of which we are interested in learning using data from experiments. This system is a common model of a servo mechanism, a mechanical system that uses feedback to control of a motor to accurately position a load mass. Examples of this include:

• Industrial robots
• 3D printers
• Tracking cameras and telescopes
• Overhead cranes

The mass that can be excited, representing the moving mass of the actuator, is equipped with a position sensor. However, the driven mass, the mass of the load we are interested in positioning, is not. In order to obtain sufficiently informative data, the second mass is instead temporarily equipped with an accelerometer. Extra temporary instrumentation is common in a lab setting, where the additional sensors are used to identify a model or to verify the performance of a state estimator etc.

A complicating factor in a setup like this is that it is difficult to accurately estimate both masses $m$ and spring stiffness $k$ at the same time. Typically, only the ratio of them is easily identifiable. This can be understood from looking at the expression for the resonance frequency $\omega$ of a simple mass-spring system, $\omega = \sqrt{k/m}$. However, if we perform two different experiments, where we change the mass by a known value, i.e., by adding a known quantity of additional mass, we can tease $k$ and $m$ apart. We will leverage this method here in order to accurately estimate all parameters of the system. We have the following parameters to estimate:

• Masses, $m_1, m_2$
• Transmission spring stiffness, $k$
• Transmission damping, $b$

The two experiments will be carried out as follows:

• Experiment 1 drives the unloaded system with a known force $F$.
• Experiment 2 adds load $\Delta m$ to the second mass, i.e., the effective mass of the load is this $m_2 + \Delta m$ where $m_2$ is unknown while $\Delta m$ is known. The same force is applied to the system as in Experiment 1.

From the classical Newton's 2nd Law of Motion, we can generate the following set of equations from the conservation equations:

\[\begin{aligned} m_1 \ddot{x_1} + b_1 \dot{x_1} &= k (x_2 - x_1) + b_2 (\dot{x_2} - \dot{x_1}) + F_1 \\ (m_2 + Δm) \ddot{x_2} + b_2 (\dot{x_2} - \dot{x_1}) + k (x_2 - x_1) &= 0.0 \\ \end{aligned}\]

Julia Environment

For this example, we will need the following packages:

using JuliaSimModelOptimizer
using ModelingToolkit
using ModelingToolkitStandardLibrary.Blocks: RealInput, TimeVaryingFunction
using OrdinaryDiffEq
using CSV, DataFrames
using DataSets
import DataInterpolations: CubicSpline
using Plots

Model Setup

Since the spring stiffness is numerically very large compared to the other parameters, and our apriori uncertainty about this parameter spans potentially several orders of magnitude, we let $log(k)$ be the parameter optimized by the solver. This brings the parameter into a similar numerical range as the other parameters, and makes the optimization problem more well-conditioned.

tvec = collect(0:0.01:200)
input_func_1 = CubicSpline(cos.(0.5 * tvec .^ 2), tvec) # A simple periodic driving force with increasing frequency

function coupledsystem(Δm = 0.0)
    @variables t x1(t) = 0.5 x2(t) = 0.5  v1(t) = 0 v2(t) = 0 mass1a(t) = 0 mass2a(t) = 0
    @parameters k = 2700.0 b1 = 9.0 b2 = 15.0 m1 = 1.5 m2 = 3.0 Δm = Δm
    D = Differential(t)
    @named F1 = RealInput()
    eqs = [
        m1 * mass1a + b1 * D(x1) ~ k * (x2 - x1) + b2 * (D(x2) - D(x1)) + F1.u,
        (m2 + Δm) * mass2a + b2 * (D(x2) - D(x1)) + k * (x2 - x1) ~ 0,
        D(x1) ~ v1,
        D(x2) ~ v2,
        D(v1) ~ mass1a,
        D(v2) ~ mass2a
    ]
    @named coupsys = ODESystem(eqs, t; systems=[F1])
    return coupsys
end

function SystemModel(; Δm = 0.0, name = :model)
    f1 = input_func_1
    coupsys = coupledsystem(Δm)
    t = ModelingToolkit.get_iv(coupsys)
    src1 = TimeVaryingFunction(f1; name=:src1)
    eqs = [ModelingToolkit.connect(coupsys.F1, src1.output)]
    ODESystem(eqs, t; systems=[coupsys, src1], name)
end

model1 = SystemModel(; Δm = 0.0)
sys1 =  complete(structural_simplify(model1))
for i in ModelingToolkit.missing_variable_defaults(sys1)
    push!(ModelingToolkit.defaults(sys1), i)
end

model2 = SystemModel(; Δm = 10.0) # We add 10kg of mass to the load in the second experiment
sys2 =  complete(structural_simplify(model2))
for i in ModelingToolkit.missing_variable_defaults(sys2)
    push!(ModelingToolkit.defaults(sys2), i)
end

Data Setup

Next, we will load some data. This data is generated from simulations. For increased realism, we simulated measurement noise on both velocity and acceleration measurements. Accelerometers have a tendency to be fairly noisy, reflected in the large amount of noise added to the acceleration measurements.

data1 = open(IO, dataset("juliasimtutorials/coupled_experiments_data_1")) do io
    CSV.read(io, DataFrame)
end

data2 = open(IO, dataset("juliasimtutorials/coupled_experiments_data_2")) do io
    CSV.read(io, DataFrame)
end

Now, lets plot the data.

plot(data1[:, "timestamp"], data1[:, "coupsys.v1"], label = "Noisy velocity measurements", layout=(2, 1), size = (1000, 600), sp=1)
plot!(data1[:, "timestamp"], data1[:, "coupsys.mass2a"], label = "Noisy acceleration measurements", sp=2)

Defining Experiments and InverseProblem

We will set up both the Experiment corresponding to Δm present or absent with their corresponding data. We will also use PredictionErrorMethod for simulations to be guided by the data which helps in the calibration process.

experiment1 = Experiment(data1, sys1, loss_func = meansquaredl2loss, alg = Tsit5(), model_transformations = [PredictionErrorMethod(0.05)])
experiment2 = Experiment(data2, sys2, loss_func = meansquaredl2loss, alg = Tsit5(), model_transformations = [PredictionErrorMethod(0.05)])

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)

Next, we specify an InverseProblem of interest to us where we specify the parameters we want to recover and their bounds. We pass both the experiments as a vector.

prob = InverseProblem([experiment1, experiment2], nothing, [sys1.coupsys.m1 => (0.0, 10.0), sys1.coupsys.m2 => (0.0, 10.0), sys1.coupsys.k => (1e1, 1e4, :log), sys1.coupsys.b1 => (0.0, 20.0), sys1.coupsys.b2 => (0.0, 20.0)])

InverseProblem with 2 experiments optimizing 5 parameters.

Calibration

Now, lets use SingleShooting for calibration. To do this, we first define an algorithm alg and then call calibrate with the prob and alg.

alg = SingleShooting(maxiters = 30, optimizer = IpoptOptimizer(tol = 1e-5, acceptable_tol = 1e-4))

r = calibrate(prob, alg)

Calibration result computed in 1 minute and 45 seconds. Final objective value: 0.00507537.
Optimization ended with Success Return Code and returned.

┌────────────┬────────────┬───────────┬────────────┬────────────┐
│ coupsys₊m1 │ coupsys₊m2 │ coupsys₊k │ coupsys₊b1 │ coupsys₊b2 │
├────────────┼────────────┼───────────┼────────────┼────────────┤
│    1.00237 │    3.99773 │    2508.0 │    10.0146 │    10.0984 │
└────────────┴────────────┴───────────┴────────────┴────────────┘

We can see the calibrated parameters closely match with their true values, i.e, k = 2500, b1 = 10.0, b2 = 10.0, m2 = 4.0.

Visualization

Lets also plot the simulation with the calibrated parameters.

plot(experiment1, prob, r, show_data = true, layout = (2, 1), size = (1000, 600), ms = 0.1)

plot(experiment2, prob, r, show_data = true, layout = (2, 1), size = (1000, 600), ms = 0.1)

We can see it fits the data well!