Custom Loss Functions

This tutorial assumes that you have read the getting started tutorial. By default JuliaSimModelOptimizer uses the squaredl2loss, but this may always not be the best choice. In this tutorial, we will work with the same model as the getting started tutorial but use a custom loss function for minimization.

Julia Environment

For this tutorial, we will use the following packages:

using JuliaSimModelOptimizer
using ModelingToolkit
using ModelingToolkit: t_nounits as t
using ModelingToolkitStandardLibrary.Electrical
using ModelingToolkitStandardLibrary.Blocks: Sine
using OrdinaryDiffEq
using Statistics
using DataSets
using Plots

Model Setup

function create_model(; C₁ = 3e-5, C₂ = 1e-6)
    @named resistor1 = Resistor(R = 5.0)
    @named resistor2 = Resistor(R = 2.0)
    @named capacitor1 = Capacitor(C = C₁)
    @named capacitor2 = Capacitor(C = C₂)
    @named source = Voltage()
    @named input_signal = Sine(frequency = 100.0)
    @named ground = Ground()
    @named ampermeter = CurrentSensor()

    eqs = [connect(input_signal.output, source.V)
        connect(source.p, capacitor1.n, capacitor2.n)
        connect(source.n, resistor1.p, resistor2.p, ground.g)
        connect(resistor1.n, capacitor1.p, ampermeter.n)
        connect(resistor2.n, capacitor2.p, ampermeter.p)]

    @named circuit_model = ODESystem(eqs, t,
        systems = [
            resistor1, resistor2, capacitor1, capacitor2,
            source, input_signal, ground, ampermeter,
        ])
end
model = create_model()
sys = complete(structural_simplify(model))

\[ \begin{align} \frac{\mathrm{d} \mathtt{capacitor2.v}\left( t \right)}{\mathrm{d}t} &= \mathtt{capacitor1.vˍt}\left( t \right) \\ 0 &= - \mathtt{resistor2.i}\left( t \right) + \mathtt{capacitor2.i}\left( t \right) - \mathtt{resistor1.i}\left( t \right) + \mathtt{capacitor1.i}\left( t \right) \end{align} \]

Data Setup

However, in this tutorial, the data is noisier. Furthermore, the strength of the noise seems to depend on the absolute value of the current $i$ going through the ampermeter. The data can be read in using the dataset function.

data = dataset("juliasimtutorials/circuit_data_heterogeneous")

name = "juliasimtutorials/circuit_data_heterogeneous"
uuid = "0a1054b8-ab9c-49d2-893a-d71525018965"
tags = []
description = "Data for the Model Optimizer loss tutorial."
type = "Blob"

[storage]
driver = "JuliaHubDataRepo"
bucket_region = "us-east-1"
bucket = "internaljuilahubdatasets"
credentials_url = "https://internal.juliahub.com/datasets/0a1054b8-ab9c-49d2-893a-d71525018965/credentials"
prefix = "datasets"
vendor = "aws"
type = "Blob"

    [storage.juliahub_credentials]
    auth_toml_path = "/home/github_actions/actions-runner-1/home/.julia/servers/internal.juliahub.com/auth.toml"

    [[storage.versions]]
    blobstore_path = "u1"
    date = "2023-08-17T01:07:47.62+00:00"
    size = 2776
    version = 1

For heteroscedastic noise, a weighted variant of l2loss outperforms the default l2loss. Specifically, we assume that the standard deviation of the measurement noise increases linearly as i increases.

Implement Custom Loss Function

We can implement such a weighted loss ourselves using the loss_func keyword of experiment. loss_func should always be a function with 3 inputs, the values of the tuned parameter or initial conditions values, the solution of the model, at the point in the search space given by the first argument and the data. Here the solution will always only contain the states and timepoints present in the dataset, e.g. state voltages will not be present in the solution. The dataset is transformed from a DataFrame to a Matrix. Both inputs are ordered in the same way, i.e. solution[i,j] corresponds to the same state and timepoint as dataset[i,j]. We can then easily implement losses using broadcasting.

function weighted_loss(tuned_vals, solution, dataset)
    σ² = var(dataset)
    sum(((solution .- dataset) .^ 2) ./ σ²)
end

weighted_loss (generic function with 1 method)

Defining Experiment and InverseProblem

In order to create an Experiment, we will use the default initial values of the states and parameters of our model. These are our initial guesses which will be used to optimize the inverse problem in order to fit the given data. We also pass in the custom loss function that we defined above.

@unpack capacitor1, capacitor2 = model
experiment = Experiment(data, sys; initial_conditions = [capacitor2.v => 0.0, capacitor1.i => 0.0], alg = Rodas4(), abstol = 1e-6, reltol = 1e-5, loss_func = weighted_loss)

Experiment for circuit_model 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, 0.1)

Once we have created the experiment, the next step is to create an InverseProblem. This inverse problem, requires us to provide the search space as a vector of pairs corresponding to the parameters that we want to recover and the assumption that we have for their respective bounds.

prob = InverseProblem(experiment, [capacitor2.C => (1.e-7, 1e-3)])

InverseProblem with one experiment with 1 elements in the search space.

Calibration

Now, lets use SingleShooting for calibration. To do this, we first define an algorithm alg and then call calibrate with the prob and alg.

alg = SingleShooting(maxiters = 10^3)
r = calibrate(prob, alg)

Calibration result computed in 18 seconds and 812 milliseconds. Final objective value: 12.0297.
Optimization ended with Success Return Code and returned.

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

Visualization

Let us visualize the calibrated results. We can plot the simulation using the calibrated parameters and compare it against the data.

plot(experiment, prob, r, show_data = true, legend = true)

We see that this calibrated current fits the data well!