Gradual parameter changes via callbacks

Specific Domains Modelling & Simulations
, ,
1

Problem

I’m trying to change the parameter of a dynamical system gradually.

In the documentation for the DifferentialEquations library there is an emphasis on abrupt / instantaneous change - when explaining callbacks ; i.e. something hits a wall; gets kicked, fast drug dosage kinetics…

However the systems I’m trying to model show gradual changes in parameters (mostly linear with time) at discrete timepoints.

While I could achieve such change using a ContinuousCallback, I’m suspecting that I should use the VectorContinuousCallback instead - however it has sadly no example yet.

Or maybe even use the DiffEqCallbacks.jl library instead?

Reproducable Example:

Gradual breaking of an harmonic oscillator, i.e. the damping force \omega gets linearly increased by 0.5 twice, at t = 25 and t = 75 each over a duration of 10 seconds.

using OrdinaryDiffEq, Plots, ModelingToolkit
# Example dynamical system:
harmonicoscillator(ddu,du,u,ω,t) =  ddu .= -ω^2 * u

# Initial Conditions, timescale, stepsize
x₀ = [0.0]; dx₀ = [π/2]; ω = .75
tspan = (0.0, 100.0); dt = 0.001

# Im using ModelingToolkit for my models mostly

@named sys = modelingtoolkitize(
  SecondOrderODEProblem(harmonicoscillator, dx₀, x₀, tspan, ω)
  )

prob = ODEProblem(complete(sys))

as previous mentioned I am looking to increasing the force ω gradually over 10 (seconds) twice;
once at t=25 and once at t=75

timeranges = [25.0:dt:35.0, 75.0:dt:85.0] # should i collect the vector here?
Δω = 0.5 # absolute change of force ω 
duration = length(timeranges[1]) # duration in timesteps of the acceleration

# This is the expected force after two changes:
expected_force = ω + ( length(timeranges) * Δω ) # i.e. :  0.75 + (2*0.5) = 1.75

# This is the naive callback i coded:

# check if the current timestep falls in one of the timeranges:
condition(u, t, integrator) =  !any(t .∈ timeranges) |> Int # hit if: == zero

# [NOTE] modelingtoolkitize assigns symbols for parameters ie force ω is now α :
affect!(integrator) = integrator.ps[:α] += Δω / duration # or (Δω * dt) / 10  
# how would i handle different durations here for example ?!

cb = ContinuousCallback(condition, affect!)

solving the system:
for the tstops argument collect the values of the ranges, vector them, splat and concatenate, I feel this is not how its supposed to go?!

sol = solve(prob, callback = cb, dt = dt, tstops = vcat(collect.(timeranges)...) ) 

plot(sol , idxs = 2); vline!(vcat(collect.(extrema.(timeranges))...))
prob.ps[:α] ≈ expected_force # true; 1.750000000000166

Question

How bad is this? how to use VectorContinuousCallback or DiffEqCallbacks.jl to achieve the same behaviour (ideally with cleaner code, and better performance).

plot output for comparison.
plot output for comparison.1002×604 73.1 KB

If this is already answered, please feel free to just point me to the resource and close this problem.
Thanks!

2

I think this is probably a bit more complicated than it needs to be. If omega only depends on time, you can just make a function that depends on time to calculate it:

harmonicoscillator(ddu,du,u,ω,t) =  ddu .= -calc_ω(t)^2 * u

function calc_ω(t)
    ω0 = 0.75
    t_dur = 10.0
    mag_increase = 0.5
    if t < 25.0
        return ω0
    elseif 25.0 ≤ t < 25.0 + t_dur 
        ω = ω0 + (t-25.0)*mag_increase/t_dur
        return ω
    elseif 25.0 + t_dur ≤ t < 75.0
        ω = ω0 + mag_increase
        return ω
    elseif 75.0 ≤ t < 75.0 + t_dur 
        ω = ω0 + mag_increase + (t-75.0)*mag_increase/t_dur
        return ω
    else 
        ω = ω0 + 2*mag_increase
        return ω
    end
end

There is probably a slicker way of writing that function, but it should at least work. You may still want to put in stops for the times (25,35,75,85) where the omega function has a discontinuous derivative.

3 Likes
4

Thanks for the answer, (i used the ifelse form, since t is a symbol not a number)

However using ifelse it works and is significantly faster than my solution.

function calc_ω(t)
    ω0 = 0.75
    t_dur = 10.0
    mag_increase = 0.5

    ω = ifelse(t < 25.0, ω0,
            ifelse(25.0 ≤ t < (25.0 + t_dur), ω0 + (t - 25.0) * mag_increase / t_dur,
            ifelse((25.0 + t_dur) ≤ t < 75.0, ω0 + mag_increase,
            ifelse(75.0 ≤ t < (75.0 + t_dur), ω0 + mag_increase + (t - 75.0) * mag_increase / t_dur,
            ω0 + 2 * mag_increase))))
    
    return ω
end

@register_symbolic calc_ω(t)

harmonicoscillator_2(ddu,du,u,ω,t) =  ddu .= -calc_ω(t)^2 * u

@named sys_2 = modelingtoolkitize(
  SecondOrderODEProblem(harmonicoscillator_2, dx₀, x₀, tspan, ω)
  )

prob_2 = ODEProblem(complete(sys_2)) 

# using callbacks;
@time sol = solve(prob, callback = cb, dt = dt, tstops = vcat(collect.(timeranges)...) ) 
#  0.125140 seconds (1.63 M allocations: 60.592 MiB)

# using the custom function
@time sol_2 = solve(prob_2, dt = dt)
#  0.000298 seconds (2.76 k allocations: 122.641 KiB)

# using the custom function with the stops:
 @time sol_2 = solve(prob_2, dt = dt, tstops = vcat(collect.(timeranges)...) )
#  0.015755 seconds (362.61 k allocations: 16.547 MiB)

So its not only due to the additional stops ..

However in this approach i wont be able to access or as a symblic parameter anymore, which makes it impossible to verify that the change worked in the intended way, and also complicates all subsequent analysis.

I still think theres a better solution using VectorContinuousCallback…

5

As you described, callbacks are for discrete changes in either state or parameters. Even VectorContinuousCallback is not applicable here, because the thing that’s continuous there is the mechanism for finding the right time to apply some change, not that change itself.

It looks like you can just make ω a LinearInterpolation from DataInterpolations.jl, which seems to have proper MTK integration: Using with Symbolics/ModelingToolkit · DataInterpolations.jl. It’s probably still a good idea to add the tstops.