Exogenous functions and data in ODE aren't working if modertly complicated

Specific Domains Numerics
,
1

Yesterday I posted about getting control inputs into an ODE and I got great advice. However, I’m running into an issue where the exogenous functions used don’t seem to work if they’re even moderately complicated.

Below is a minimal example. I have an ODE defined that pulls in a constant c and exogenous function D_(t). If D(t) uses a simple condition, e.g., time > 5, it works as intended. If I add anything slightly more complicated, like another condition, e.g., (time > 5) && (time < 6), it doesn’t appear to have any effect anymore. In my real project I’m trying to pull in values in D from data, but even just hard-coding a control signal does not seem to work. I’m not sure where I’m going wrong in this.

Minimal working example:

using DifferentialEquations
using ParameterizedFunctions
using Optim
using Plots
gr()

function D1(t)
    if (t > 5)
        return 0.1
    else
        return 0.0
    end
end
function D2(t)
    if (t > 5) && (t < 6)
        return 0.1
    else
        return 0.0
    end
end
p = [
    -0.1089 ;
    75.2304 ;
     2.3965
]
c = p[end]

sma_ode1 = @ode_def begin
    dT = a*T + b*D1(t) + c
end a b
sma_ode2 = @ode_def begin
    dT = a*T + b*D2(t) + c
end a b
u0 = [22.0]
p0 = p[1:2]
tspan = (0.0, 10.0)
prob1 = ODEProblem(sma_ode1, u0, tspan, p0)
prob2 = ODEProblem(sma_ode2, u0, tspan, p0)
sol1 = solve(prob1, Tsit5(), saveat=0.1)
sol2 = solve(prob2, Tsit5(), saveat=0.1)
p11 = plot(sol1, ylims=(0,50), title="Single condition")
p12 = plot(sol1.t, D1.(sol1.t))
p21 = plot(sol2, ylims=(0,50), title="Two conditions")
p22 = plot(sol2.t, D2.(sol2.t))
l = @layout [a b ; c d]
plot(p11, p21, p12, p22, layout = l)

Output:

imagen
imagen719×486 9.42 KB

1 Like
2

The adaptive-timestep ODE methods take larger timesteps in smooth regions. In the first case, the forcing function stays on, and the solver notices that the forcing function has changed and backtracks to better capture the change. In the second case, the forcing function returns to zero, and the timestep was large enough to step right past the active time of the forcing function. You can force the solver to stop at the discontinuous timesteps with tstops=[5.0, 6.0], decrease the maximum timestep, or change the forcing function with a callback.

3

Awesome! Thank you very much for your help! So in the case of data (which looks similar, as the control signal for this experiment is largely just switching on and off) I suppose what I should do is set tstop to be points of discontinuity I detect just using the derivative of the control.

4

Note that saveat doesn’t force the ODE solver to step through all those times - it just interpolates the solution:
image

With tstops=[5, 6],
image

5

If you don’t mind a related followup: I’m actually trying to use DiffEqParamEstim to estimate parameters a and b. Are you aware if this time step behavior affects optimizations as well? I don’t see a similar tstops parameter from the optimize routine or problem setup. Maybe this needs to become a separate question.

6

I haven’t used DiffEqParamEstim, but looking at the docs, you should be able to pass additional keyword arguments for the ODE solver when constructing the loss function.

7

And Handling Instability When Solving ODE Problems - #5 by ChrisRackauckas goes into some strategies for handling discontinuous inputs.

1 Like