3

Similar to this question, I am trying to solve this ODE with a time-dependent input parameter. It consists of a series of discrete callbacks. At certain times, a parameter is changed (not a state!). Times and values are stored in a nx2 Array. But I can't get the affect function to find the corresponding parameter value at the specified time. In the given examples, the value assigned to u[1] is usually constant. Consider this MWE (with a very Matlab-like approach), which works correctly without the callback:

using DifferentialEquations
using Plots

function odm2prod(dx, x, params, t)
    k_1, f_1, V_liq, X_in, Y_in, q_in = params

    rho_1 = k_1*x[1]
    q_prod = 0.52*f_1*x[1]
    # Differential Equations
    dx[1] = q_in/V_liq*(X_in - x[1]) - rho_1
    dx[2] = q_in/V_liq*(Y_in - x[2])
end

x0      = [3.15, 1.5]
tspan   = (0.0, 7.0)
params  = [0.22, 43, 155, 249, 58, 0]
prob    = ODEProblem(odm2prod, x0, tspan, params)

input   = [1.0 60; 1.1 0; 2.0 60; 2.3 0; 4.0 430; 4.05 0]
dosetimes = input[:,1]
function affect!(integrator)
    ind_t = findall(integrator.t == dosetimes)
    integrator.p[6] = input[ind_t, 2]
end
cb = PresetTimeCallback(dosetimes, affect!)
sol = solve(prob, Tsit5(), callback=cb, saveat=1/12)

plot(sol, vars=[1, 2])

It does not work. The error originates at line 22, since comparing a vector to a scalar seems not to be defined in Julia, or there is a special syntax I am unaware of.

I know that it is possible to use time-dependent parameters in Julia, but I suppose that would only work for continuous functions, not discrete changes!? I haven taken a look at the help for interpolate, but I am not sure how to use it for my specific case.

Could someone tell me how to get this to work, please? Should probably need just a few lines of code. Also, I do not necessarily want dosetimes as part of sol.t, unless they coincide.

Chris Rackauckas
  • 18,645
  • 3
  • 50
  • 81
winkmal
  • 622
  • 1
  • 6
  • 16

1 Answers1

3

You are using findall wrong, the documentation says

findall(f::Function, A)

Return a vector I of the indices or keys of A where f(A[I]) returns true.

Then you have to take into account that the result of a search for "all" is a list. As you expect it to only have one element, use the first one only

function affect!(integrator)
    ind_t = findall(t -> t==integrator.t, dosetimes)
    integrator.p[6] = input[ind_t[1], 2]
end

and you get the plot

plot of the ode solution

Community
  • 1
  • 1
Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
  • Actually, I wanted to mimic Matlab's `ismember` function with Julia's `findall`. But then, I should have placed the arguments vice-versa... Anyway, your solutions works. An alternative would be to use Interpolations.jl, but the necessary option for `previous` interpolation is [still pending](https://github.com/JuliaMath/Interpolations.jl/pull/338). – winkmal Feb 21 '20 at 20:58
  • I was surprised to see this work even when the first `PresetTimeCallback` happens at the initial time point, e.g. `tspan[1]`. I am working on a problem where `PresetTimeCallback` modifies the state and not a parameter. In that case, the callback does not get triggered at the initial time point. Any idea why this would make a difference? – fabern Jul 02 '20 at 15:05
  • 1
    Anyway, if it helps someone: in the 2nd case where `affect!` modifies a state, e.g. `integrator.u[1] += input[ind_t[1], 2]` a workaround could consist of two steps: a) use `PresetTimeCallback(dosetimes, affect!, initialize = (c,u,t,integrator) -> affect!(integrator))` (as mentioned in https://github.com/SciML/DifferentialEquations.jl/issues/570#issuecomment-605385113). and b) modify `affect!` in order to check `if (integrator.t in dosetimes)`, as otherwise `findall` might throw an error. – fabern Jul 02 '20 at 15:08
  • Hi @fabern, didn't have a look here for some time now. I was coding my first steps in Julia, to see if I could change to it permanently. But this jump feature is an essential aspect in the kind of modeling I do, so I need a really solid solution for it. It seems that [Pumas](https://docs.pumas.ai/dev/tutorials/introduction/#Building-a-Subject-1) would be an interesting alternative, since their `DosageRegimen` very much resembles the kind of feed I use. I will give it a try eventually. – winkmal Aug 03 '20 at 13:13