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:

ModuleDescription
JuliaSimModelOptimizerThis is used to formulate our inverse problem and solve it
FMI and FMIImportLibrary for importing and simulating FMUs
CSV and DataFramesWe will read our experimental data from .csv files
PlotsThe 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
RowtimestampJ1.phiJ1.phi.derJ2.phiJ2.phi.derJ3.phiJ3.phi.derJ4.phiJ4.phi.der
Float64Float64Float64Float64Float64Float64Float64Float64Float64
10.00.010.00.00.00.00.00.00.0
20.010.09955279.915840.0004993590.0997420.00.00.00.0
30.020.1984219.862870.001989630.1979210.00.00.00.0
40.030.2969059.838220.00444740.2929890.00.00.00.0
50.040.3952629.83650.007834020.3834460.00.00.00.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 48 seconds.
┌────────┬──────────┐
│ freqHz │       T2 │
├────────┼──────────┤
│    2.0 │ 0.999999 │
│    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 │      1.0 │
│    2.0 │      1.0 │
│    2.0 │ 0.999999 │
│    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))
Example block output

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.996439 │
│ 2.00371 │ 0.995232 │
│ 2.00371 │  0.99641 │
│ 2.00371 │ 0.995134 │
│ 2.00371 │ 0.998241 │
│ 2.00371 │ 0.994298 │
│ 2.00371 │ 0.994727 │
│ 2.00371 │ 0.999141 │
│ 2.00371 │  0.99683 │
│ 2.00371 │ 0.994347 │
│ 2.00371 │  0.99396 │
│ 2.00371 │ 0.993809 │
│ 2.00371 │ 0.990773 │
│ 2.00371 │ 0.994095 │
│ 2.00371 │ 0.992065 │
│ 2.00371 │ 0.993779 │
│    ⋮    │    ⋮     │
└─────────┴──────────┘
        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))
Example block output

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