Consider the following example:
import Plots
using ModelingToolkit, DifferentialEquations
@parameters t σ ρ β
@variables x(t) y(t) z(t)
@derivatives D'~t
eqs = [D(x) ~ σ * (y - x),
D(y) ~ x * (ρ - z) - y,
D(z) ~ x * y - β * z]
lorenz1 = ODESystem(eqs, name=:lorenz1)
lorenz2 = ODESystem(eqs, name=:lorenz2)
@variables a
@parameters γ
connections = [0 ~ lorenz1.x + lorenz2.y + a * γ]
connected = ODESystem(connections, t, [a], [γ], systems=[lorenz1,lorenz2])
u0 = [lorenz1.x => 1.0,
lorenz1.y => 0.0,
lorenz1.z => 0.0,
lorenz2.x => 0.0,
lorenz2.y => 1.0,
lorenz2.z => 0.0,
a => 2.0]
The parameters are than specified as a vector of pairs, corresponding to the parameters of the ODESystem.
p = [lorenz1.σ => 10.0,
lorenz1.ρ => 28.0,
lorenz1.β => 8 / 3,
lorenz2.σ => 10.0,
lorenz2.ρ => 28.0,
lorenz2.β => 8 / 3,
γ => 2.0]
tspan = (0.0, 100.0)
prob = ODEProblem(connected, u0, tspan, p)
sol = solve(prob, Rodas5())
It feels a bit inconsistent to speficy p as a vector of pairs but inside the integrator it is just a vector.
In a callback one would have to use e.g. int.p[1] which in this case corresponds to gamma.
julia> prob.p
7-element Array{Float64,1}:
2.0
10.0
28.0
2.6666666666666665
10.0
28.0
2.6666666666666665
It seems to be constructed like this:
[
connected.ps...,
connected.systems[1].ps...,
connected.systems[2].ps...,
]
Is there an easier way of referencing the parameters e.g. for usage in a callback function?