Calibrate a Functional Mock-up Unit (FMU)
In this example, we will calibrate a Functional Mock-up Unit (FMU), following Functional Mockup Interface v2.0, to experiment data using the SingleShooting
method.
The example is based on a Coupled-Clutches model that is available as an FMU. This Coupled-Clutches FMU is a model of a clutch system with three clutches.
Julia environment
For this example, we will use the following packages:
Module | Description |
---|---|
JuliaSimModelOptimizer | This is used to formulate our inverse problem and solve it |
FMI and FMIImport | Library for importing and simulating FMUs |
CSV and DataFrames | We will read our experimental data from .csv files |
Optimization | High level interface package for using optimizers |
OptimizationNLopt | The optimization solver package with GN_ISRES |
Plots | The plotting and visualization library |
using JuliaSimModelOptimizer
using FMI, FMIImport
using DataFrames, CSV
using Optimization, OptimizationNLopt
using Plots
Model Setup
We start off by loading the FMU using FMI.jl - an open-source julia package which enables working with FMUs. We specify the type of the FMU as FMI.fmi2TypeModelExchange
since we intend to simualate the FMU in Model Exchange mode.
fmu = fmiLoad(joinpath(@__DIR__, "../assets/CoupledClutches_ME.fmu");
type = FMI.fmi2TypeModelExchange)
Model name: CoupledClutches
Type: 0
Note that since the FMU support needs both FMI and FMIImport, both packages need to be loaded.
Data Setup
We load the data from a CSV file using the DataFrames.jl package. The data is a time series of the states of the system, such as the angular velocity and acceleration of the inertias. The true values of the parameters are freqHz = 2.0
and T2 = 1.0
.
data = CSV.read(joinpath(@__DIR__, "../assets/CoupledClutches_data.csv"), DataFrame)
Row | timestamp | J1.phi | J1.phi.der | J2.phi | J2.phi.der | J3.phi | J3.phi.der | J4.phi | J4.phi.der |
---|---|---|---|---|---|---|---|---|---|
Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | |
1 | 0.0 | 0.0 | 10.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
2 | 0.01 | 0.0995527 | 9.91584 | 0.000499359 | 0.099742 | 0.0 | 0.0 | 0.0 | 0.0 |
3 | 0.02 | 0.198421 | 9.86287 | 0.00198963 | 0.197921 | 0.0 | 0.0 | 0.0 | 0.0 |
4 | 0.03 | 0.296905 | 9.83822 | 0.0044474 | 0.292989 | 0.0 | 0.0 | 0.0 | 0.0 |
5 | 0.04 | 0.395262 | 9.8365 | 0.00783402 | 0.383446 | 0.0 | 0.0 | 0.0 | 0.0 |
Defining Experiment and InverseProblem
We next use the data to setup an experiment. We specify the time span of the Experiment
and the relative tolerance of the solver.
experiment = Experiment(data, fmu; tspan = (0.0, 1.5), reltol = 1e-8)
Experiment for Functional Mock-up Unit 2.0 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, 1.5)
Next, We define an InverseProblem
of interest to us where we specify the parameters we want to recover and their bounds.
prob = InverseProblem(experiment, fmu, ["freqHz" => (0.05, 3.0), "T2" => (0.2, 1.5)])
InverseProblem with one experiment optimizing 2 parameters.
Calibration
We run the calibration using the calibrate
function. We specify the SingleShooting
algorithm to use for the calibration. We also specify the adtype
to be NoAD()
as we will use a gradient free optimizer for the calibration process. We use GN_ISRES
as the optimizer.
opt = OptimizationNLopt.NLopt.GN_ISRES()
alg = SingleShooting(maxiters = 30_000, optimizer = opt)
r = calibrate(prob, alg, adtype = Optimization.SciMLBase.NoAD())
r
Calibration result computed in 3 minutes and 48 seconds. Final objective value: 1.31002e-30.
Optimization ended with MaxIters Return Code and returned.
┌────────┬─────┐
│ freqHz │ T2 │
├────────┼─────┤
│ 2.0 │ 1.0 │
└────────┴─────┘
We can see the values of the calibration matches the true values of the parameters, i.e., freqHz = 2.0
and T2 = 1.0
from the data.
Visualization
We can now plot the simulation using the calibrated parameters and compare it against the data.
plot(experiment, prob, r, show_data = true, legend = :best, ms = 0.2, layout = (8, 1), size = (1000, 1600), left_margin = 5Plots.mm)