Using optimization parameters in ModelingToolkit

Specific Domains Modelling & Simulations
1

Hi All,

I want to be able to pass optimization parameters to the model.

I took the RC circuit from the docs and created 2 parameters (param_r1 param_c1)
https://mtk.sciml.ai/stable/tutorials/acausal_components/#Connecting-and-Simulating-Our-Electrical-Circuit

My code:

using ModelingToolkit, DifferentialEquations
using ModelingToolkitStandardLibrary.Electrical

@variables t
@parameters param_r1 param_c1

@named resistor = Resistor(R=param_r1)
@named capacitor = Capacitor(C=param_c1)
@named source = ConstantVoltage(V=1.)
@named ground = Ground()

rc_eqs = [
    connect(source.p, resistor.p)
    connect(resistor.n, capacitor.p)
    connect(capacitor.n, source.n)
    connect(capacitor.n, ground.g)
]

@named _rc_model = ODESystem(rc_eqs, t)
@named rc_model = compose(_rc_model,
    [resistor, capacitor, source, ground])
sys = structural_simplify(rc_model)
u0 = [
    capacitor.v => 0.0
]

params = [param_r1 => 1.0, param_c1 => 1.0]
tspan = (0., 10.)

prob = ODAEProblem(sys, u0, tspan, params)

The last line gives me the next error:

ERROR: MethodError: no method matching float(::Type{Any})
Closest candidates are:
  float(::Type{Union{Missing, T}}) where T at C:\Users\mzhen\AppData\Local\Programs\Julia-1.7.2\share\julia\base\missing.jl:112       
  float(::Any) at C:\Users\mzhen\AppData\Local\Programs\Julia-1.7.2\share\julia\base\float.jl:269
  float(::Union{StatsBase.PValue, StatsBase.TestStat}) at C:\Users\mzhen\.julia\packages\StatsBase\pJqvO\src\statmodels.jl:86
  ...
Stacktrace:
 [1] float(#unused#::Type{Any})
   @ Base .\missing.jl:113
 [2] promote_to_concrete(vs::Vector{Any}; tofloat::Bool, use_union::Bool)
   @ ModelingToolkit C:\Users\mzhen\.julia\packages\ModelingToolkit\LXLjs\src\utils.jl:558
 [3] varmap_to_vars(varmap::Vector{Pair{Num, Float64}}, varlist::Vector{Sym{Real, Base.ImmutableDict{DataType, Any}}}; defaults::Dict{Any, Any}, check::Bool, toterm::Function, promotetoconcrete::Nothing, tofloat::Bool, use_union::Bool)
   @ ModelingToolkit C:\Users\mzhen\.julia\packages\ModelingToolkit\LXLjs\src\variables.jl:68
 [4] ODAEProblem{true}(sys::ODESystem, u0map::Vector{Pair{Num, Float64}}, tspan::Tuple{Float64, Float64}, parammap::Vector{Pair{Num, Float64}}; callback::Nothing, use_union::Bool, check::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})     
   @ ModelingToolkit.StructuralTransformations C:\Users\mzhen\.julia\packages\ModelingToolkit\LXLjs\src\structural_transformation\codegen.jl:523
 [5] ODAEProblem{true}(sys::ODESystem, u0map::Vector{Pair{Num, Float64}}, tspan::Tuple{Float64, Float64}, parammap::Vector{Pair{Num, Float64}})
   @ ModelingToolkit.StructuralTransformations C:\Users\mzhen\.julia\packages\ModelingToolkit\LXLjs\src\structural_transformation\codegen.jl:514
 [6] ODAEProblem(::ODESystem, ::Vararg{Any}; kw::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ ModelingToolkit.StructuralTransformations C:\Users\mzhen\.julia\packages\ModelingToolkit\LXLjs\src\structural_transformation\codegen.jl:492
 [7] ODAEProblem(::ODESystem, ::Vararg{Any})
   @ ModelingToolkit.StructuralTransformations C:\Users\mzhen\.julia\packages\ModelingToolkit\LXLjs\src\structural_transformation\codegen.jl:492
 [8] top-level scope
   @ d:\Projects\playground\julia\grid\src\grid _stl.jl:31

How to fix it?

Thanks!

2

@YingboMa this looks like a bug?

3

Don’t use new parameters like param_r1 and param_c1 while resistor.R and capacitor.R already parameters. MTK doesn’t know param_r1 and param_c1 are parameters of the model.

4

Thank you @YingboMa.

I fixed my code:

using ModelingToolkit, Plots, DifferentialEquations
using ModelingToolkitStandardLibrary.Electrical

@variables t

@named resistor = Resistor(R=1.)
@named capacitor = Capacitor(C=1.)
@named source = ConstantVoltage(V=1.)
@named ground = Ground()

rc_eqs = [
    connect(source.p, resistor.p)
    connect(resistor.n, capacitor.p)
    connect(capacitor.n, source.n)
    connect(capacitor.n, ground.g)
]

@named _rc_model = ODESystem(rc_eqs, t)
@named rc_model = compose(_rc_model,
    [resistor, capacitor, source, ground])
sys = structural_simplify(rc_model)
u0 = [
    capacitor.v => 0.0
]

params = [resistor.R => 1.0, capacitor.C => 1.0]
params_init = [1.0, 1.0]
tspan = (0., 10.)

prob = ODAEProblem(sys, u0, tspan, params)
sol = solve(prob, Tsit5(), p=params_init)

But the solution is wrong (constant zeros)

julia> sol
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 9-element Vector{Float64}:
  0.0
  1.0e-6
  1.1e-5
  0.00011099999999999999
  0.0011109999999999998
  0.011110999999999996
  0.11111099999999996
  1.1111109999999995
 10.0
u: 9-element Vector{Vector{Float64}}:
 [0.0]
 [1.3014e-320]
 [1.4313e-319]
 [1.44425e-318]
 [1.4448273e-317]
 [1.4377578e-316]
 [1.36834641e-315]
 [8.72596552e-315]
 [-1.78158667517e-312]
5

@YingboMa Any idea why the solution is wrong? ^^^^

6

It should be

sol = solve(prob, Tsit5())
7

@YingboMa
But I need to run optimization in respect of parameters p. It means I need to call the function solve with the different parameters many times. It seems solve supports such API.
For example in the docs, here you can see the call:

solve(prob_sde, SOSRI(), p=p,
               sensealg = ForwardDiffSensitivity(), saveat = 0.1)

When I do

params = [resistor.R => 2.0, capacitor.C => 2.0]
prob = ODAEProblem(sys, u0, tspan, params)

julia> prob.p
3-element Vector{Float64}:
 2.0
 2.0
 1.0

So instead of 2 parameters, we have 3, with remake we have 2 parameters.

prob1 = remake(prob; p=[2., 2.])
julia> prob1.p
2-element Vector{Float64}:
 2.0
 2.0

But the solution is still wrong:

sol = solve(prob1, Tsit5())
julia> sol
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 7-element Vector{Float64}:
  0.0
  9.999999999999999e-5
  0.0010999999999999998
  0.011099999999999997
  0.11109999999999996
  1.1110999999999995
 10.0
u: 7-element Vector{Vector{Float64}}:
 [0.0]
 [0.0]
 [0.0]
 [0.0]
 [0.0]
 [0.0]
 [0.0]

julia>
8

For now (until we fix it), you need to be more careful with remake on MTK generated problems. See the FAQ about this: Frequently Asked Questions · ModelingToolkit.jl

9

@ChrisRackauckas @YingboMa

It seems my problem is in the set of parameters to set up.
If I use all 3 parameters to remake and solve, everything is OK.
But if I want to play only with 2 of them, it gives me the wrong solution.

So question is - what is the right way to play with only some subset of parameters and keep other defaults?

For example:
sol = solve(prob1, Tsit5(), p=[1., 3., 7])
gives the correct result

sol = solve(prob1, Tsit5(), p=[1., 3.])
gives the wrong result.

10

Can you open an issue? It’s weird it doesn’t just error with the wrong number of parameters so it would be good to dig into this one in more depth.

11

@ChrisRackauckas
I’ve opened the issue.
https://github.com/SciML/ModelingToolkit.jl/issues/1631

But what is the right way of solving the system with the different parameters (a subset of parametrs) supposed to be?

1 Like
12

@ChrisRackauckas @YingboMa
Sorry to bother you guys, but I just want to know the “official way” to run Solver with the subset of parameters, is it possible?

If it’s not supported, what is the right way to do optimization? Every time update only “optimization parameters” and keep the rest unmodified?

Thanks!

13

Your case has 3 parameters. Passing 2 just gives undefined behavior. We’re using this as a test case to add a size check in the MTK generated functions and in the remake phase to throw earlier. But in general, the Pair constructors are the only ones you really should be using here since otherwise you could be making assumptions about parameter ordering which may not be true post-simplification.