Skip to content

Calibration Analysis

For many models, we don't know all the parameters when creating them, so we need to calibrate them and find the optimal parameters. To do that, we can either rely on experimental data that we have for the modelled system, or we can tune the parameters such that the system achieves a certain state. In this tutorial, we will fit a model such that the result of the simulation matches some experimental data.

Method Overview

The Dyad compiler uses the physical description of the system to generate a system of differential-algebraic equations:

x=f(x,y,p,t)0=g(x,y,p,t)

If we call z(t,p) to be the subset of the state [x(t,p),y(t,p)] for which we have data d(ti) at time points ti, the goal is to find the parameters p such that z(ti,p)d(ti). This is done by defining a loss function L(p) such that p is the argument which minimizes the L. A nonlinear optimization is then performed to find the minimizer of L. In particular, gradient-based optimization using automatic differentiation and adjoint methods is used to do this fast and robustly.

Example Definition

To run a calibration analysis, we need a model and some data that corresponds to the behavior that we would want from it. In this example we will solve the inverse problem of finding the model parameters with the CalibrationAnalysis, starting from an electrical circuit called a capacitor bridge, also known as a de Sauty bridge.

Adding required depencies

This example requires the following julia packages to be installed in the current environment. Use the Julia REPL to add them:

julia
using Pkg
Pkg.add(["ElectricalComponents", "BlockComponents", "DyadData", "DyadModelOptimizer"])

Let's first create the component to analyze, in Dyad. If you already have a component library, you can add this into a new .dyad file in the dyad folder, otherwise you need to create a new component library via the Command Palette in VS Code (⇧ ⌘ P on MacOS, Ctrl-Shift-P on Windows). See the getting started instructions for more details.

Dyad code

We can now add the following Dyad code

dyad
example component MyDeSauty
  resistor1 = ElectricalComponents.Resistor(R=5.0)
  resistor2 = ElectricalComponents.Resistor(R=2.0)
  capacitor1 = ElectricalComponents.Capacitor(C=C₁)
  capacitor2 = ElectricalComponents.Capacitor(C=C₂)
  current_sensor = ElectricalComponents.CurrentSensor()
  source = ElectricalComponents.VoltageSource()
  input_signal = BlockComponents.Sine(frequency=100.0, amplitude=1, phase=0)
  ground = ElectricalComponents.Ground()
  parameter C₁::Capacitance = 30μ
  parameter C₂::Capacitance = 45μ
relations
  initial capacitor2.v = 0.0
  connect(input_signal.y, source.V)
  connect(source.p, capacitor1.n, capacitor2.n)
  connect(source.n, resistor1.p, resistor2.p)
  connect(source.n, ground.g)
  connect(resistor1.n, current_sensor.n, capacitor1.p)
  connect(resistor2.n, current_sensor.p, capacitor2.p)
end

The next thing that we need is a dataset. In this case we will use an example dataset that is hosted on JuliaHub via DyadData.jl. With DyadData.jl we can represent a timeseries-like dataset that is backed by a JuliaHub dataset via the DyadDataset type.

DyadData.DyadDataset Method
julia
DyadDataset(
username::AbstractString,
name::AbstractString;
independent_var::AbstractString,
dependent_vars::Vector{<:AbstractString},
kwargs...)

Represent a timeseries-like dataset that is backed by a JuliaHub dataset.

Keyword arguments

  • independent_var: the name of the column that represents the independent variable (usually the time)

  • dependent_vars: a vector of the names of the columns for the dependent variables

When reading files (local file option or a downloaded JuliaHub dataset), CSV.jl is used. Additional keyword arguments passed to this function will be passed on to CSV.read. This can help with changing settings such as the delimiter used in the file. See https://csv.juliadata.org/stable/reading.html for more details.

source

Now we can write a calibration analysis for this MyDeSauty component, extending the CalibrationAnalysis from DyadModelOptimizer.jl. Since this functionality is provided by DyadModelOptimizer.jl, you first have to explicitly add that package and also DyadData for specifying datasets:

dyad
using DyadModelOptimizer: CalibrationAnalysis
using DyadData: DyadDataset

Next we are pass all the keywords of the calibration analysis to the CalibrationAnalysis:

dyad
analysis DeSautyCalibrationAnalysis
  extends DyadModelOptimizer.CalibrationAnalysis(
    data=DyadData.DyadDataset("juliasimtutorials", "circuit_data", independent_var = "timestamp", dependent_vars = ["ampermeter.i(t)"]),
    N_cols=1, depvars_cols=["ampermeter.i(t)"], depvars_names=["current_sensor.i"],
    N_tunables = 1, search_space_names=["capacitor2.C"],
    search_space_lb=[1e-7], search_space_ub=[1e-3],
    calibration_alg="SingleShooting", optimizer_maxiters=100)
  model = MyDeSauty()
end

Let's compare the behaviour of the system before an after calibration. To check the original system behaviour, we can make use of the TransientAnalysis to simulate the system in time.

dyad
analysis DeSautyTransient
  extends TransientAnalysis(abstol=1e-8, reltol=1e-8, stop=0.1)
  model = MyDeSauty()
end

Now, let's simulate and plot!

julia
using DyadInterface, DyadData, Plots

spec = DeSautyTransientSpec()
sol0 = run_analysis(spec)
data = DyadDataset("juliasimtutorials", "circuit_data", independent_var = "timestamp", dependent_vars = ["ampermeter.i(t)"])
df = build_dataframe(data)
plot(sol0, idxs=spec.model.current_sensor.i)
scatter!(df.timestamp, df."ampermeter.i(t)")

Finally, we can run the calibration analysis:

julia
res = DeSautyCalibrationAnalysis()
Calibration Analysis solution for CalibrationAnalysis: Calibration result computed in 22 seconds and 419 milliseconds. Final objective value: 2.63095e-5.
Optimization ended with Success Return Code and returned.

┌──────────────┐
│ capacitor2₊C │
├──────────────┤
│   1.00065e-5 │
└──────────────┘

and then plot the calibration results. In this case, we retrieve an artifact that contains the plot.

julia
artifacts(res, :CalibrationResultPlotExperiment1)

We can also view the numeric values of the parameters via another artifact:

julia
artifacts(res, :CalibrationResultTableExperiment1)
1×1 DataFrame
Rowcapacitor2₊C
Float64
11.00065e-5

Analysis Arguments

Required Arguments

  • name: The name of the analysis.

  • model: An ODESystem representing the model that will be used for numerical integration.

  • stop: The end time for the integration.

  • data: A DyadDataset from DyadData.

  • indepvar: The independent variable in the data that was passed.

  • N_cols: The number of dependent variable columns in the data.

  • depvars_cols: A vector of column names for dependent variables.

  • depvars_names: A vector of names in the model that map to the depvars_cols. Note that they need to be in the same order.

  • N_tunables: The number of tunable parameters.

  • search_space_names: The names of the tunable parameters (as a vector of strings).

  • search_space_lb: A vector of values for the lower bound of the tunable parameters.

  • search_space_ub: A vector of values for the upper bound of the tunable parameters.

  • optimizer_maxiters: The maximum number of iterations for the optimizer.

  • optimizer_maxtime: The maximum number of (real-world) seconds for the optimizer to run.

Optional Arguments

  • alg: The ODE integrator to use as a symbol. Possible options are: "auto" (default), "Rodas5P", "FBDF", "Tsit5".

  • start: The start time for the integration. Defaults to 0.0.

  • abstol: Absolute tolerance to use during the simulation. Defaults to 1e-8.

  • reltol: Relative tolerance to use during the simulation. Defaults to 1e-8.

  • saveat: Timepoints to save the solution at or 0 (to let the integrator decide, default).

  • dtmax: The maximum allowed timestep or 0 (to let the integrator decide, default).

  • loss_func: The loss function to use for comparing data to simulation results. Defaults to "l2loss". Options are: "l2loss", "norm_meansquaredl2loss", "meansquaredl2loss, "squaredl2loss".

  • calibration_alg: The calibration algorithm to use. Available options are: "SingleShooting", "MultipleShooting", "SplineCollocation" & "KernelCollocation".

  • multiple_shooting_trajectories: If MultipleShooting is selected, specify the number of trajectories.

  • pem_gain: Gain factor for Prediction Error Method. If this is greater than 0, the DiscreteFixedGainPEM model transformation is applied.

  • optimizer: The optimizer to use. Available options are: "auto", "Ipopt", "BBO".

  • optimizer_abstol: The absolute tolerance to use for the optimizer. Defaults to 1e-4.

Artifacts

A CalibrationAnalysis returns the following artifacts:

Customizable Plot

The customizable plot of a CalibrationAnalysis is the timeseries plot corresponding to the simulation of the system with the parameters given by the results of the calibration compared against the data.

Standard Plots

  • CalibrationResultPlotExperiment1: Returns the standard time series plot of the solution trajectory for the calibrated parameters overlayed with the given data points.

  • ConvergencePlot: Returns a plot of the loss function evolution during calibration with a logarithmic y axis.

Tables

  • CalibrationResultTableExperiment1: Returns a DataFrame of the calibrated parameters.

Further Reading

Layout Switch

Adjust the layout style of VitePress to adapt to different reading needs and screens.

Expand all
The sidebar and content area occupy the entire width of the screen.
Expand sidebar with adjustable values
Expand sidebar width and add a new slider for user to choose and customize their desired width of the maximum width of sidebar can go, but the content area width will remain the same.
Expand all with adjustable values
Expand sidebar width and add a new slider for user to choose and customize their desired width of the maximum width of sidebar can go, but the content area width will remain the same.
Original width
The original layout width of VitePress

Page Layout Max Width

Adjust the exact value of the page width of VitePress layout to adapt to different reading needs and screens.

Adjust the maximum width of the page layout
A ranged slider for user to choose and customize their desired width of the maximum width of the page layout can go.

Content Layout Max Width

Adjust the exact value of the document content width of VitePress layout to adapt to different reading needs and screens.

Adjust the maximum width of the content layout
A ranged slider for user to choose and customize their desired width of the maximum width of the content layout can go.

Spotlight

Highlight the line where the mouse is currently hovering in the content to optimize for users who may have reading and focusing difficulties.

ONOn
Turn on Spotlight.
OFFOff
Turn off Spotlight.