ModelingToolkitStandardLibrary and Optimization

Specific Domains Modelling & Simulations
1

Hi All,

I want to optimize the Model created based on ModelingToolkitStandardLibrary.

So I have a small RC-Model copied from the official documentation to ModelingToolkitStandardLibrary.

module Grid

export get_model

using ModelingToolkit
using ModelingToolkitStandardLibrary.Electrical
using ModelingToolkitStandardLibrary.Blocks: Constant

function get_model(; resistance::Float64=1.0, capacitance::Float64=1.0, voltage::Float64=1.0)
    @variables t
    @named R = Resistor(R=resistance)
    @named C = Capacitor(C=capacitance)
    @named V = Voltage()
    @named ground = Ground()
    @named source = Constant(k=voltage)

    rc_eqs = [
        connect(V.p, R.p)
        connect(R.n, C.p)
        connect(C.n, V.n, ground.g)
        connect(V.V, source.output)
    ]

    @named _rc_model = ODESystem(rc_eqs, t)

    @named rc_model = compose(_rc_model, [R, C, V, source, ground])
    sys = structural_simplify(rc_model)
    return sys, states(rc_model)
end

end

The full code to run optimization is:

using Revise

includet("./grid.jl")
using .Grid

using ModelingToolkit, DifferentialEquations
using Optimization, OptimizationOptimisers, DiffEqSensitivity, Zygote, OptimizationOptimJL
using CSV, Plots, DataFrames, Statistics

function get_obs(fname::String, obs_vars::Vector{String})
    df = DataFrame(CSV.File(fname))
    return df[:, "time"], df[:, obs_vars]
end

function get_param_vars(model::ModelingToolkit.AbstractSystem)
    params = parameters(model)
    return params
end

function loss(p::Vector{Float64})
    pred_sol = solve(prob, Tsit5(), p=p, saveat=t_save)
    pred = Array(pred_sol')
    mse = mean((pred - obs_vals) .^ 2)
    return mse
end

# observation variable names
obs_vars = ["C₊v(t)"]
obs_data_file = "data/res.csv"

t_save, obs = get_obs(obs_data_file, obs_vars)
model, model_state = get_model()
param_vars = get_param_vars(model)

u0 = [s => 0.0 for s in states(model)]
tspan = (t_save[1], last(t_save))
lb = [0.0, 0.0, 0.99]
ub = [3.0, 3.0, 1.01]

prob = ODAEProblem(model, u0, tspan)
obs_vals = Array(obs)

p_init = [1.0, 1.0, 1.0]

adtype = Optimization.AutoZygote()
optf = Optimization.OptimizationFunction((x, p) -> loss(x), adtype)
optprob = Optimization.OptimizationProblem(optf, p_init, lb=lb, ub=ub)
result = Optimization.solve(optprob, LBFGS(), maxiters=300)

That code works well, and I can get the correct solution.
My problem is in the ordering of the parameters. The official doc says “symbol ordering is not guaranteed” but if I do something like this:

function loss(p::Vector{Float64})

    p1 = [param_vars[1] => p[1], param_vars[2] => p[2], param_vars[3] => p[3]]
    p2 = ModelingToolkit.varmap_to_vars(p1, param_vars)

    pred_sol = solve(prob, Tsit5(), p=p2, saveat=t_save)
    pred = Array(pred_sol')
    mse = mean((pred - obs_vals) .^ 2)
    return mse
end

I get the next error:

ERROR: Compiling Tuple{Type{Dict}, Vector{Pair{Sym{Real, Base.ImmutableDict{DataType, Any}}, Float64}}}: try/catch is not supported.
Refer to the Zygote documentation for fixes.

So how can I do optimization and guarantee the correct order of the parameters?

2

Can you open an issue? We probably need to add an adjoint to varmap_to_vars to avoid the try/catch

3

I’ve created the issue:
https://github.com/SciML/ModelingToolkit.jl/issues/1692

4

@ChrisRackauckas
Two more questions:

  1. Any workaround until it gets fixed?
  2. There are three parameters: R₊R, C₊C, and source₊k, but we know the source, so we don’t need to optimize it. How we can play just with two parameters (R₊R, C₊C) instead of three.

Thanks!

5

I posted a workaround in the issue. See if that works and if it does then we should add it to the FAQ.

Make loss(p) with a vector of 2 things and put the constant inside of the loss function.

6

Thank you @ChrisRackauckas

I responded here:

I want to understand the source of uncertainty, If I run params = parameters(model) and get
params = [param_1, param_2, ..., param_n]

And then I do solve(prob, Tsit5(), p=[p_1, p_2, ..., p_n])

Does it mean that it is possible to have: p_i != param_i?

7

yes