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:
If we call
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.
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 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.
artifacts(res, :CalibrationResultPlotExperiment1)
We can also view the numeric values of the parameters via another artifact:
artifacts(res, :CalibrationResultTableExperiment1)
Row | capacitor2₊C |
---|---|
Float64 | |
1 | 1.00065e-5 |
Analysis Arguments
Required Arguments
name
: The name of the analysis.model
: AnODESystem
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 thedepvars_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 to1e-8
.reltol
: Relative tolerance to use during the simulation. Defaults to1e-8
.saveat
: Timepoints to save the solution at or0
(to let the integrator decide, default).dtmax
: The maximum allowed timestep or0
(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
: IfMultipleShooting
is selected, specify the number of trajectories.pem_gain
: Gain factor for Prediction Error Method. If this is greater than 0, theDiscreteFixedGainPEM
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 to1e-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.