# 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:

Module | Description |
---|---|

JuliaSimModelOptimizer | This is used to formulate our inverse problem and solve it |

ModelingToolkit | The symbolic modeling environment |

ModelingToolkitStandardLibrary | Library for using standard modeling components |

OrdinaryDiffEq | The numerical differential equation solvers |

Statistics | Library for standard statistics functions |

DataSets | We will load our experimental data from datasets on JuliaHub |

Plots | The plotting and visualization library |

```
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!