# Using optimization parameters in ModelingToolkit

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!

@YingboMa this looks like a bug?

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.

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

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

It should be

``````sol = solve(prob, Tsit5())
``````

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

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

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.

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.

@ChrisRackauckas
Iâ€™ve opened the issue.

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

1 Like

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

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.