I may not be thinking about this right, but currently I have a system defined as a heirarchy of models. I store the full system in the parameters of a DifferentialEquations.jl ODEProblem. I can then use methods directly on the system or models to calculate the ode function. The system is variable in its size and configuration at the start of the sim, but does not change once the sim is running. I initialize the ODEProblem by looking into the system struct and defining p and u0 based on the characteristics of the system.
An example of one of the models could be a step input signal. Normally I would implement a step input as DiscreteCallback, working through the integrator interface. I would like to set the field for the step input model in the parameters, always referencing the same field for a given callback.
Is there a way to somehow define a callback so that it always references some field in the parameters, without explicitly calling out the parameter itself?
Attempt at a MWE below.
using DifferentialEquations
mutable struct Model
initial_state::Float64
signal_value::Float64
end
function odefunc!(du,u,p,t)
for i in eachindex(p.models)
du[i] = p.models[i].signal_value * t
end
return nothing
end
function initialize_models(models)
u = getfield.(models,:initial_state)
p = (;models = [model1,model2])
for model in models
# define callbacks for each model
# in models that will change the field
# signal_value of model in the parameters p
# to be used in the odefunc calculation
# put the callback into a CallbackSet to be used in solve
# affect!(integrator) = integrator.p.( somehow always reference model1 ).signal_value = 1.0
end
return u,p#,cbs
end
model1 = Model(1.0, 0.0)
model2 = Model(0.0, 1.0)
models = [model1,model2]
u0,p = initialize_models(models)
#u0,p,cbs = initialize_models(models)
prob = ODEProblem(odefunc!, u0, (0.,10.), p)
sol = solve(prob) #cbs
I’m not entirely sure what you’re trying to do, but did you try using an anonymous function affect! = (integrator) -> integrator.p.model[models].signal_value = 1.0? Defining dispatches of functions in a loop isn’t a great idea, while an anonymous function will act more naturally.
1 Like
Thanks @ChrisRackauckas I had thought about using anonymous functions but I guess I wasn’t sure the right way to do it at the time. Looking at your code, I would assume [models] is globally scoped or I guess function scoped. I wasn’t sure if that would work or if it was the right way to do this. I think in the end I just over thought it, and it looks like your suggestion will work for me.
An actual working MWE is below that maybe better indicates what I’m trying to do. It seems to be working when I’m fixed step but doesn’t seem to be producing the right results when variable step (obviously a separate issue).
using DifferentialEquations
mutable struct Model
initial_state::Float64
step_time::Float64
signal_value::Float64
Model(initial_state,step_time) = new(initial_state,step_time,0.0)
end
function odefunc!(du,u,p,t)
for i in eachindex(p.models)
du[i] = p.models[i].signal_value * 1.0 #arbitrary function
end
return nothing
end
function initialize_models(models)
u = getfield.(models,:initial_state)
p = (;models = models)
cbs = []
for i in eachindex(models)
push!(cbs, DiscreteCallback(
(u,t,integrator) -> t >= integrator.p.models[i].step_time, #condition
(integrator) -> integrator.p.models[i].signal_value = 1.0) #affect!
)
end
cbs = CallbackSet(cbs...)
return u,p,cbs
end
function simulate(original_models, adaptive = true)
# copy since we're going to mutate models during sim
# this way we can can rerun simulate(models)
models = deepcopy(original_models)
u0,p,cbs = initialize_models(models)
prob = ODEProblem(odefunc!, u0, (0.,10.), p)
if adaptive
sol = solve(prob, Tsit5(), callback = cbs)
else
sol = solve(prob, Tsit5(), adaptive = false, dt = 1, callback = cbs)
end
return sol
end
model1 = Model(0.0,1.0)
model2 = Model(-1.0,2.0)
models = [model1,model2]
#sol = simulate(models, true) #doesnt seem to work
sol = simulate(models, false) #works
variable step:
fixed step:
Note that discrete callbacks apply at the end of steps. t >= integrator.p.models[i].step_time means it will apply at the first step which ends up past integrator.p.models[i].step_time. It will not apply at integrator.p.models[i].step_time. If you want that to happen, then you should set a tstop (or see PresetTimeCallback)
1 Like
Yup, that’ll do it. Missed that detail. Thanks again Chris.
