Parametric Uncertainty Quantification with Functional Mockup Units (FMUs)

This example demonstrates how use Functional Mockup Units (FMUs), following Functional Mockup Interface v2.0, to perform parametric uncertainty quantification with JuliaSimModelOptimizer.

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
Plots	The plotting and visualization library

using JuliaSimModelOptimizer
using FMI, FMIImport
using CSV, DataFrames
using Plots

Model Exchange

Model Setup

Let us first demonstrate the use of the Coupled-Clutches FMU in Model Exchange mode.

We start of 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

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.

data = CSV.read(joinpath(@__DIR__, "../assets/CoupledClutches_data.csv"), DataFrame)

5×9 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 in the Experiment and the relative tolerance of the solver.

experiment = Experiment(data, fmu; tspan = (0.0, 1.5), reltol = 1e-5)

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 specify an InverseProblem of interest to us where we specify the parameters we want to recover and their bounds.

invprob = InverseProblem(experiment, fmu, ["freqHz" => (0.05, 3), "T2" => (0.2, 1.5)])

InverseProblem with one experiment optimizing 2 parameters.

Parametric Uncertainty Quantification

We now run the parameteric uncertainty quantification using the parametric_uq function. We specify the InverseProblem and the StochGlobalOpt algorithm. We also specify the sample size for the analysis.

results = parametric_uq(invprob, StochGlobalOpt(; maxiters = 3000), sample_size = 20)

Parametric uncertainty ensemble of length 20 computed in 7 minutes and 58 seconds.
┌────────┬──────────┐
│ freqHz │       T2 │
├────────┼──────────┤
│    2.0 │      1.0 │
│    2.0 │      1.0 │
│    2.0 │      1.0 │
│    2.0 │      1.0 │
│    2.0 │      1.0 │
│    2.0 │      1.0 │
│    2.0 │      1.0 │
│    2.0 │      1.0 │
│    2.0 │      1.0 │
│    2.0 │      1.0 │
│    2.0 │      1.0 │
│    2.0 │ 0.999999 │
│    2.0 │      1.0 │
│    2.0 │      1.0 │
│    2.0 │      1.0 │
│    2.0 │      1.0 │
│   ⋮    │    ⋮     │
└────────┴──────────┘
       4 rows omitted

Note

Note that we need a sample size of at least 5000 to get a good result and this low sample size is just for demonstration. Same with the number of maxiters parameter of StochGlobalOpt.

Visualization

We can see the results of the analysis by plotting its histogram.

plot(results, layout = (2, 1), bins = 8, show_comparison_value = false, legend = :best, size = (1000, 1000))

We can see the distribution for both the parameters is concentrated near their true values.

Co-Simulation Mode

Model Setup

Let us now import the FMU in Co-Simulation mode.

fmu = fmiLoad(joinpath(@__DIR__, "../assets/CoupledClutches_CS.fmu"); type = FMI.fmi2TypeCoSimulation)

Model name:        CoupledClutches
Type:              1

Defining Experiment and InverseProblem

The steps for using FMUs in Co-Simulation mode are similar to those for Model Exchange mode. The only difference is that we specify the type of the FMU as FMI.fmi2TypeCoSimulation and we use the DiscreteExperiment instead of the Experiment to setup the experiment as Co-Simulation FMUs are discrete time models and specify the step size (here 1e-2) compulsorily.

de = DiscreteExperiment(data, fmu, 1e-2; tspan = (0.0, 1.5))
invprob = InverseProblem(de, fmu, ["freqHz" => (0.05, 3), "T2" => (0.2, 1.5)])

InverseProblem with one experiment optimizing 2 parameters.

Parametric Uncertainty Quantification

The rest of the code remains the same as for the Model Exchange FMU.

results = parametric_uq(invprob, StochGlobalOpt(; maxiters = 3000), sample_size = 20)

Parametric uncertainty ensemble of length 20 computed in 1 minute and 57 seconds.
┌─────────┬──────────┐
│  freqHz │       T2 │
├─────────┼──────────┤
│ 2.00371 │ 0.991535 │
│ 2.00371 │ 0.993547 │
│ 2.00371 │ 0.995275 │
│ 2.00371 │  0.99424 │
│ 2.00371 │ 0.992794 │
│ 2.00371 │ 0.994766 │
│ 2.00371 │ 0.990951 │
│ 2.00371 │ 0.991749 │
│ 2.00371 │ 0.998018 │
│ 2.00371 │ 0.993346 │
│ 2.00371 │  0.99118 │
│ 2.00371 │ 0.994102 │
│ 2.00371 │ 0.996396 │
│ 2.00371 │ 0.995778 │
│ 2.00371 │ 0.996437 │
│ 2.00371 │ 0.991155 │
│    ⋮    │    ⋮     │
└─────────┴──────────┘
        4 rows omitted

Visualization

We can see the results of the analysis by plotting its histogram.

plot(results, layout = (2, 1), bins = 8, show_comparison_value = false, legend = :best, size = (1000, 1000))

Again, we can see the distribution for both the parameters is concentrated near their true values.