Getting Started with JuliaSimModelOptimizer

Accessing JuliaSimModelOptimizer

In order to start using JuliaSimModelOptimizer, you can either sign in on JuliaHub and start the JuliaSim IDE application or use the JuliaHubRegistry to download the package locally.

Using JuliaSimModelOptimizer on JuliaHub

After you sign in on JuliaHub (or your JuliaHub enterprise instance), go to Applications and then start a JuliaSim IDE app.

Using JuliaSimModelOptimizer locally

In order to use the JuliaSimModelOptimizer locally, you'll need to first add PkgAuthentication to your global environment.

using Pkg
original_active_project = Base.active_project()
try
    withenv("JULIA_PKG_SERVER" => nothing) do
        Pkg.activate("v$(VERSION.major).$(VERSION.minor)"; shared = true)
        Pkg.add("PkgAuthentication")
    end
finally
    Pkg.activate(original_active_project)
end

Next, configure Julia to be aware of JuliaHub's package server (or your enterprise JuliaHub instance):

using PkgAuthentication
PkgAuthentication.install("juliahub.com")

Then add all known registries on JuliaHub, including the JuliaHub registry, from the "Pkg REPL" (by pressing ] from the Julia REPL), using

Pkg.Registry.add()

This will open up a new browser window or tab and ask you to sign in to JuliaHub or to confirm authenticating as an already logged in user.

Lastly, in order to use the datasets in the tutorials, you'll have to make them accessible via JuliaHubData, so you'll have to install JuliaHubData and DataSets and then run the following:

using JuliaHubData
using DataSets
using PkgAuthentication

PkgAuthentication.authenticate()
JuliaHubData.connect_juliahub_data_project()

Fitting model parameters to data with JuliaSimModelOptimizer

For many of the models that we have, we don't know all the parameters, so we need to calibrate them. 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. We will first look at the case where we have some experimental data and we want to fit the model such that the simulation matches the data. In the following tutorial we will solve the inverse problem of finding the model parameters starting from a capacitor bridge, also sometimes called De Sauty bridge.

This capacitor bridge consists of 2 capacitors, 2 resistors and an AC source with a frequency of 100Hz and 1V amplitude. We assume that we know the resistor values and the value of one of the capacitances, as shown above. We want to find the value of the unknown capacitance starting from some measurements obtained from the ampermeter.

First steps

For this tutorial we will need the following packages:

using JuliaSimModelOptimizer
using ModelingToolkit, OrdinaryDiffEq
using ModelingToolkitStandardLibrary.Electrical
using ModelingToolkitStandardLibrary.Blocks: Sine
using DataSets
using Plots

Model Setup

In order to solve the inverse problem, we must first define the circuit model.

function create_model(; C₁ = 3e-5, C₂ = 1e-6)
    @variables t
    @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 = structural_simplify(model)

\[ \begin{align} \frac{\mathrm{d} capacitor2_{+}v\left( t \right)}{\mathrm{d}t} =& capacitor1_{+}vˍt\left( t \right) \\ 0 =& - resistor2_{+}i\left( t \right) + capacitor2_{+}i\left( t \right) - resistor1_{+}i\left( t \right) + capacitor1_{+}i\left( t \right) \end{align} \]

In the above, we first define the circuit components using ModelingToolkitStandardLibrary and then we connect them as in the diagram above. The units for the model are all in SI, as indicated in its documentation. The default value that we have provided for $C_2$ is our initial guess for the capacitance value.

Data Setup

We can read tabular experimental data where the model names are matching column names in the table. The time dependence in the column names is optional. For example let's consider the following:

data = dataset("juliasimtutorials/circuit_data")

name = "juliasimtutorials/circuit_data"
uuid = "56ac293a-9ebf-401c-b8a9-57d986fb2747"
tags = []
description = "Data for the Model Optimizer getting started tutorial."
type = "Blob"

[storage]
driver = "JuliaHubDataRepo"
bucket_region = "us-east-1"
bucket = "internaljuilahubdatasets"
credentials_url = "https://internal.juliahub.com/datasets/56ac293a-9ebf-401c-b8a9-57d986fb2747/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-06-10T23:53:27.845+00:00"
    size = 2768
    version = 1

We can now create an Experiment that would correspond to the experimental data above. The inverse problem is defined by the experiment, the model of the system and the search space. The search space is a vector of pairs, where the first element is the model variable (parameter or initial condition) that we want to find and the second represents the bounds in which we believe the value can be.

@unpack capacitor2 = model
experiment = Experiment(data, sys; alg = Rodas4(), abstol = 1e-6, reltol = 1e-5)
prob = InverseProblem(experiment, sys, [capacitor2.C => (1.e-7, 1e-3)])

InverseProblem with one experiment optimizing 1 parameters.

Calibration

We can now calibrate the model parameters to the data using a calibration method of our choice. The simplest method that we can start with is SingleShooting. We use IpoptOptimizer for defining our optimizer.

alg = SingleShooting(maxiters = 100, optimizer = IpoptOptimizer(; tol = 1e-6))
r = calibrate(prob, alg)

Calibration result computed in 19 seconds and 319 milliseconds. Final objective value: 7.09128e-10.
Optimization ended with Success Return Code and returned.

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

To check that the result is correct, we can plot the solution of the simulated ampermeter readings to the experimental ones. To include the data in the plot, we use show_data=true. As systems can have a very large amount of states, the legend is not shown by default, so we turn that on with legend=true.

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

We can also try another calibration algorithm, such as MultipleShooting, which would split the data into multiple groups and fit them individually and impose a continuity condition via a penalty in the loss function, for which the magnitude is set with continuitylossweight. We can also specify maxtime as stopping criteria instead of maxiters.

alg = MultipleShooting(maxiters = 200, trajectories = 5, continuitylossweight = 1.0, optimizer = IpoptOptimizer(; tol = 1e-6), initialization = DataInitialization())
r = calibrate(prob, alg)

Calibration result computed in 1 minute and 54 seconds. Final objective value: 7.41535e-10.
Optimization ended with MaxIters Return Code and returned.

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

Again, to check that the solution is correct, we can plot the result of the simulated experiment with the calibration result and compare it against the data.

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

If we want to inspect other states of the model, we can use the states keyword argument.

@unpack capacitor2, resistor1 = model
plot(experiment, prob, r, states = [capacitor2.v, resistor1.v], legend = true, ylabel = "V")

Design Optimization with JuliaSimModelOptimizer

Similarly to how we can use a Wheatstone bridge to determine an unknown resistance by varying a variable resistance, we can build a capacitance bridge in which we find and unknown capacitance. In the following, instead of assuming that we can vary one of the resistances, we directly calibrate the capacitance value of the unknown capacitor, such that the current through the ampermeter is 0.

Design problems usually require tuning model parameters such that certain goals are attained. In some cases our objective function might not be depend on measurement data, such as in fitting problems, but instead we might want the system to achieve a certain state. JuliaSimModelOptimizer can be used for both of these usecases and in the following we will present a design optimization problem starting from the above defined capacitor bridge.

t0 = 0.0
tend = 0.1

@unpack ampermeter, capacitor2 = model

config = DesignConfiguration(sys; tspan = (t0, tend), alg = Rodas4(),
                        running_cost = abs2(ampermeter.i),
                        reduction = sum)

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

InverseProblem with one experiment optimizing 1 parameters.

As we don't have data in this case, we need to provide a timespan for our experiment and an error function to be minimized. The inverse problem definition looks the same as above, as we have the same system and the same search space.

Calibration

We define the experiment without any data, but we supply a custom error function, such that we enforce that no current passes through the ampermenter. This is equivalent with balancing the bridge by modyfing the resistances. This means that we can verify our calibration result with the known analytic results for this particular configuration,

\[C_2 = C_1 * \frac{R_1}{R_2}.\]

The simplest calibration method that we can try is SingleShooting.

alg = SingleShooting(maxiters = 100, optimizer = IpoptOptimizer(; tol = 1e-6))
r = calibrate(prob, alg)

Calibration result computed in 20 seconds and 467 milliseconds. Final objective value: 4.17521e-14.
Optimization ended with Success Return Code and returned.

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

As we can see, we obtain the correct result $C_2 = 75μF$. We can check that the system behaves how we expect, namely that we don't have current through the ampermeter. To do this, we can plot result of the experiment using the calibration result that we obtained above.

plot(config, prob, r, states = ampermeter.i, legend = true)

We can also try MultipleShooting, so that we split the above oscillations in multiple groups and instead of trying to cancel all of them at the same time, we can just find a solution that minimizes parts of a single oscillation and enforce the continuity of the solution.

alg = MultipleShooting(maxtime = 240, trajectories = 50, continuitylossweight = 1e-4, initialization = RandomInitialization())
r = calibrate(prob, alg)

Calibration result computed in 33 seconds. Final objective value: 0.00247788.
Optimization ended with Success Return Code and returned.

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

We can again check that the system behaves as expected.

plot(config, prob, r, states = ampermeter.i, legend = true)