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
Let us load the required packages.
using FMI
using DataFrames, CSV
using JuliaSimModelOptimizer
using Plots
Model Exchange Mode
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
Load the data
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)
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 |
Experiment Setup
We next use the data to setup an experiment. We specify the time span of the Experiment
and the relative tolerance of the solver. We also specify the InverseProblem
of interest to us where we specify the parameters we want to recover and their bounds.
ex = Experiment(data, fmu; tspan = (0.0, 1.5), reltol = 1e-5)
invprob = InverseProblem(ex, fmu, ["freqHz" => (0.05, 3), "T2" => (0.2, 1.5)])
InverseProblem for Functional Mock-up Unit 2.0 with one experiment optimizing 2 parameters.
Run 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 55 seconds.
┌────────┬──────────┐
│ freqHz │ T2 │
├────────┼──────────┤
│ 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 │
│ 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 │
│ ⋮ │ ⋮ │
└────────┴──────────┘
4 rows omitted
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
.
Visualize the results
We can see the results of the analysis by plotting its histogram.
plot(results, layout=(2,1), bins = 8, show_comparision_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
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.
fmu = fmiLoad(joinpath(@__DIR__, "../assets/CoupledClutches_CS.fmu"); type = FMI.fmi2TypeCoSimulation)
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 for Functional Mock-up Unit 2.0 with one experiment optimizing 2 parameters.
Run 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 2 minutes and 24 seconds.
┌─────────┬──────────┐
│ freqHz │ T2 │
├─────────┼──────────┤
│ 2.00371 │ 0.997265 │
│ 2.00371 │ 0.990419 │
│ 2.00371 │ 0.992624 │
│ 2.00371 │ 0.996264 │
│ 2.00371 │ 0.996497 │
│ 2.00371 │ 0.997826 │
│ 2.00371 │ 0.994068 │
│ 2.00371 │ 0.997419 │
│ 2.00371 │ 0.99467 │
│ 2.00371 │ 0.992008 │
│ 2.00371 │ 0.996819 │
│ 2.00371 │ 0.992827 │
│ 2.00371 │ 0.991196 │
│ 2.00371 │ 0.99167 │
│ 2.00371 │ 0.992533 │
│ 2.00371 │ 0.99772 │
│ ⋮ │ ⋮ │
└─────────┴──────────┘
4 rows omitted
Visualize the results
We can see the results of the analysis by plotting its histogram.
plot(results, layout=(2,1), bins = 8, show_comparision_value = false, legend = :best, size = (1000, 1000))
Again, we can see the distribution for both the parameters is concentrated near their true values.