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:

$$\begin{array}{rl}{x}^{\prime }& =f(x,y,p,t)\\ 0& =g(x,y,p,t)\end{array}$$

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(t_i)$ at time points $t_i$, the goal is to find the parameters $p\ast$ such that $z(t_i,p\ast )\approx d(t_i)$. This is done by defining a loss function $L(p)$ such that $p\ast$ 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:

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

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

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:

using DyadModelOptimizer: CalibrationAnalysis
using DyadData: DyadDataset

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

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.

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

Now, let's simulate and plot!

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:

res = DeSautyCalibrationAnalysis()

Calibration Analysis solution for CalibrationAnalysis: Calibration result computed in 31 seconds. 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.

artifacts(res, :CalibrationResultPlotExperiment1)

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

artifacts(res, :CalibrationResultTableExperiment1)

1×1 DataFrame

Row	capacitor2₊C
	Float64
1	1.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

• DyadModelOptimizer's Documentation

• Generalized physics-informed learning through language-wide differentiable programming

• Differentiable Programming for Differential Equations: A Review

• Solving Inverse Problems in Julia: Introduction to Julia's SciML from Scratch

• SciMLSensitivity.jl's Documentation