# ModelingToolkitStandardLibrary and Optimization

**URL:** <https://discourse.julialang.org/t/modelingtoolkitstandardlibrary-and-optimization/84146>\
**Category:** Modelling & Simulations\
**Tags:** modelingtoolkit\
**Created:** [July 13, 2022, 6:08am UTC](https://discourse.julialang.org/t/modelingtoolkitstandardlibrary-and-optimization/84146 "2022-07-13T06:08:40Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![mzhenirovskyy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mzhenirovskyy/32/7549_2.png) [@mzhenirovskyy](https://discourse.julialang.org/u/mzhenirovskyy)\
**Post date:** [July 13, 2022, 6:08am UTC](https://discourse.julialang.org/t/modelingtoolkitstandardlibrary-and-optimization/84146/1 "2022-07-13T06:08:41Z")

</div>

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.

```julia
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:

```julia
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:

```julia
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:

```julia
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?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 13, 2022, 8:41am UTC](https://discourse.julialang.org/t/modelingtoolkitstandardlibrary-and-optimization/84146/2 "2022-07-13T08:41:52Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![mzhenirovskyy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mzhenirovskyy/32/7549_2.png) [@mzhenirovskyy](https://discourse.julialang.org/u/mzhenirovskyy)\
**Post date:** [July 13, 2022, 7:20pm UTC](https://discourse.julialang.org/t/modelingtoolkitstandardlibrary-and-optimization/84146/3 "2022-07-13T19:20:06Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![mzhenirovskyy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mzhenirovskyy/32/7549_2.png) [@mzhenirovskyy](https://discourse.julialang.org/u/mzhenirovskyy)\
**Post date:** [July 14, 2022, 5:01am UTC](https://discourse.julialang.org/t/modelingtoolkitstandardlibrary-and-optimization/84146/4 "2022-07-14T05:01:09Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 14, 2022, 7:20am UTC](https://discourse.julialang.org/t/modelingtoolkitstandardlibrary-and-optimization/84146/5 "2022-07-14T07:20:04Z")

</div>

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

> [@mzhenirovskyy](#):
>
> 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.

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

---

<div class="post-metadata">

**Author:** ![mzhenirovskyy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mzhenirovskyy/32/7549_2.png) [@mzhenirovskyy](https://discourse.julialang.org/u/mzhenirovskyy)\
**Post date:** [July 16, 2022, 2:53am UTC](https://discourse.julialang.org/t/modelingtoolkitstandardlibrary-and-optimization/84146/6 "2022-07-16T02:53:27Z")

</div>

Thank you @ChrisRackauckas

I responded [here](https://github.com/SciML/ModelingToolkit.jl/issues/1692#issuecomment-1184797598):

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 16, 2022, 5:15am UTC](https://discourse.julialang.org/t/modelingtoolkitstandardlibrary-and-optimization/84146/7 "2022-07-16T05:15:06Z")

</div>

> [@mzhenirovskyy](#):
>
> Does it mean that it is possible to have: `p_i != param_i`?

yes