Skip to content

Commit

Permalink
add phase offset option to PeriodicCallback
Browse files Browse the repository at this point in the history
  • Loading branch information
baggepinnen committed May 28, 2024
1 parent a35955b commit b52a59f
Show file tree
Hide file tree
Showing 2 changed files with 51 additions and 11 deletions.
28 changes: 17 additions & 11 deletions src/iterative_and_periodic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -82,8 +82,6 @@ struct PeriodicCallbackAffect{A, dT, Ref1, Ref2}
end

function (S::PeriodicCallbackAffect)(integrator)
@unpack affect!, Δt, t0, index = S

add_next_tstop!(integrator, S)

affect!(integrator)
Expand All @@ -93,7 +91,7 @@ function add_next_tstop!(integrator, S)
@unpack Δt, t0, index = S

# Schedule next call to `f` using `add_tstops!`, but be careful not to keep integrating forever
tnew = t0[] + (index[] + 1) * Δt
@show tnew = t0[] + (index[] + 1) * Δt
tstops = integrator.opts.tstops
#=
Okay yeah, this is nasty
Expand All @@ -110,34 +108,42 @@ end

"""
```julia
PeriodicCallback(f, Δt::Number; initial_affect = false,
PeriodicCallback(f, Δt::Number; phase = 0, initial_affect = false,
final_affect = false,
kwargs...)
```
`PeriodicCallback` can be used when a function should be called periodically in terms of
integration time (as opposed to wall time), i.e. at `t = tspan[1]`, `t = tspan[1] + Δt`,
`t = tspan[1] + 2Δt`, and so on. This callback can, for example, be used to model a
`t = tspan[1] + 2Δt`, and so on.
If a non-zero `phase` is provided, the invokations of the callback will be shifted by `phase` time units, i.e., the calls will occur at

Check warning on line 120 in src/iterative_and_periodic.jl

View workflow job for this annotation

GitHub Actions / Spell Check with Typos

"invokations" should be "invocations".
`t = tspan[1] + phase`, `t = tspan[1] + phase + Δt`,
`t = tspan[1] + phase + 2Δt`, and so on.
This callback can, for example, be used to model a
discrete-time controller for a continuous-time system, running at a fixed rate.
## Arguments
- `f` the `affect!(integrator)` function to be called periodically
- `Δt` is the period
## Keyword Arguments
## Keyword Arguments
- `phase` is a phase offset
- `initial_affect` is whether to apply the affect at `t=0`, which defaults to `false`
- `final_affect` is whether to apply the affect at the final time, which defaults to `false`
- `kwargs` are keyword arguments accepted by the `DiscreteCallback` constructor.
"""
function PeriodicCallback(f, Δt::Number;
phase = 0,
initial_affect = false,
final_affect = false,
initialize = (cb, u, t, integrator) -> u_modified!(integrator,
initial_affect),
kwargs...)

phase < 0 && throw(ArgumentError("phase offset must be non-negative"))
# Value of `t` at which `f` should be called next:
t0 = Ref(typemax(Δt))
index = Ref(0)
Expand All @@ -154,8 +160,8 @@ function PeriodicCallback(f, Δt::Number;
initialize_periodic = function (c, u, t, integrator)
@assert integrator.tdir == sign(Δt)
initialize(c, u, t, integrator)
t0[] = t
index[] = 0
t0[] = t + phase
index[] = iszero(phase) ? 0 : -1
if initial_affect
affect!(integrator)
else
Expand Down
34 changes: 34 additions & 0 deletions test/periodic_tests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -126,3 +126,37 @@ sol2 = solve(prob, Tsit5(), callback = cb)
@test sol2(72.0 + eps(72.0)) [10.0]
@test sol2(96.0 + eps(96.0)) [10.0]
@test sol2(120.0 + eps(120.0)) [10.0]

# Test phase offset
approxin(needle, haystack) = any(isapprox(needle), haystack)

function integ(du, u, p, t)
du[1] = 1
du[2] = 0
end

function nullaffect!(integ)
integ.u[2] += 1
end

cb = PeriodicCallback(nullaffect!, 1.0; phase = 0.1)
prob = ODEProblem(integ, [0.0, 0], (0.0, 10.0))
sol = solve(prob, Tsit5(), callback = cb)
@test sol.u[end][2] == 10 # Test expected number of calls to affect
for t in 0.1:1:10
@test approxin(t, sol.t) # Test that the integrator stopped at all expected points
end
sol.t[end] 10 # test that we did not step over end

# With negative phase
@test_throws ArgumentError PeriodicCallback(nullaffect!, 1.0; phase = -0.1)

# Phase offset larger than period
cb = PeriodicCallback(nullaffect!, 1.0; phase = 1.1)
prob = ODEProblem(integ, [0.0, 0], (0.0, 10.0))
sol = solve(prob, Tsit5(), callback = cb)
@test sol.u[end][2] == 9 # Test expected number of calls to affect
for t in 1.1:1:10
@test approxin(t, sol.t) # Test that the integrator stopped at all expected points
end
sol.t[end] 10 # test that we did not step over end

0 comments on commit b52a59f

Please sign in to comment.