Callbacks for parameters in a stiff system of ODEs?

Specific Domains Modelling & Simulations
1

I have a stiff system of 21 ODEs, for which I want to use callback functions on some time-dependent parameters.

Looking for some help on this, as I’m running into some error messages and I’m new to callbacks. Not entirely sure how to incorporate the piecewise initial conditions.

Basically, I’m looking for:

v1 = 0 and v2 = 0 if t < a + d
v1 = gh and v2 = g(1-h) if t >= a+d
ydot = zeros(21,1) if t < a
ydot = ydot if t >= a

This is my code:

using LinearAlgebra, DifferentialEquations, MATLABDiffEq

a = 25;
d = 5;
T = 2323;
V = 1;
c = 0.1;
d = 0.35;
e = 0.5;
f = 0.2;
g = 1.0/3.0;
h = 0.3;
k = 0.12;
l = 1.0/11;
m = 0.6;
n = 0.1;
p = n*l;
q = (1-n)*l;
z = 0.01;
s = 0.2;
t = q;
u = 0.5*l;
v = 0.592;
w = v*u;
x = (1-v)*u;
y = 0.23245;
z  = 13;

function RHS(ydot,y,par,t)

ydot = 0.0*zeros(21,1);
y = 0.0*zeros(21,1);
par = 0.0*zeros(23,1);
v1 = 0.0;
v2 = 0.0;
a = par[1];
d = par[2];
T = par[3];
d = par[4];
c = par[5];
e = par[6];
f = par[7];
g = par[8];
h = par[9];
k = par[10];
l = par[11];
m = par[12];
n = par[13];
p = par[14];
q = par[15];
z = par[16];
s = par[17];
t = par[18];
u = par[19];
w = par[20];
x = par[21];
y = par[22];
z = par[23];

fM = d * y[1] * (y[19] + e * (y[5] + y[8] + y[11] + y[14]) + f * y[17])/T;
h = c * d * y[2] * (y[19] + e * (y[5] + y[8] + y[11] + y[14]) + f * y[17])/T;

ydot[1] = -fM - (v1*y[1]-v2*y[2]);
ydot[2] = -h + (v1*y[1]-v2*y[2]);
ydot[3] = fM - k*y[3]-(v1*y[3]-v2*y[4]);
ydot[4] = h - k*y[4]+(v1*y[3]-v2*y[4]);
ydot[5] = k*(y[3]-y[5])-z*y[5]-(v1*y[5]-v2*y[6]);
ydot[6] = k*(y[4]-y[6])-z*y[6]+(v1*y[5]-v2*y[6]);
ydot[7] = z*(y[5]+y[6])-k*y[7];
ydot[8] = k*(y[5]-y[8])-z*y[8]-(v1*y[8]-v2*y[9]);
ydot[9] = k*(y[6]-y[9])-z*y[9]+(v1*y[8]-v2*y[9]);
ydot[10] = z*(y[8]+y[9])-k*y[10] + k*y[7];
ydot[11] = k*(y[8]-y[11])-z*y[11]-(v1*y[11]-v2*y[12]);
ydot[12] = k*(y[9]-y[12])-z*y[12]+(v1*y[11]-v2*y[12]);
ydot[13] = z*(y[11]+y[12])-k*y[13] + k*y[10];
ydot[14] = k*(y[11]-y[14])-z*y[14]-(v1*y[14]-v2*y[15]);
ydot[15] = k*(y[12]-y[15])-z*y[15]+(v1*y[14]-v2*y[15]);
ydot[16] = z*(y[14]+y[15])-k*y[16] + k*y[13];
ydot[17] = m*k*y[14]-z*y[17]-(v1*y[17]-v2*y[18])-t*y[17];
ydot[18] = m*k*y[15]-z*y[18]+(v1*y[17]-v2*y[18])-t*y[18];
ydot[19] = (1-m)*k*y[14]-(z+s)*y[19]-(v1*y[19]-v2*y[20])-(p+q)*y[19];
ydot[20] = (1-m)*k*y[15]-(z+s)*y[20]+(v1*y[19]-v2*y[20])-(p+q)*y[20];
ydot[21] = (1-m)*k*(y[14]+y[15]);

return ydot
end

function condition(out,y,t,integrator) # Event when event_f(u,t) == 0
out[1] = integrator.t - integrator.par[1] - integrator.par[2];
out[2] = integrator.par[1] - integrator.t;
end

function affect!(integrator, idx)
if idx == 1
    v1 = integrator.par[8]*integrator.par[9];
    v2 = integrator.par[8]*(1-integrator.par[9]);
else
    ydot = 0.0*zeros(21,1); 
end
end

cb = VectorContinuousCallback(condition,affect!,2);
y0 = [2323 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1];
tspan = (0.0, 100.0);
prob = ODEProblem(RHS,y0,tspan);
sol = solve(prob, MATLABDiffEq.ode15s(), callback = cb, saveat = 1.0)

and the error with stack trace:

ERROR: LoadError: BoundsError
Stacktrace:
 [1] getindex at ./number.jl:78 [inlined]
 [2] RHS(::Array{Any,2}, ::Array{ModelingToolkit.Operation,2}, ::Array{Any,1}, ::ModelingToolkit.Operation) at /Users/dc/desktop/is/post.jl:60
 [3] (::ODEFunction{true,typeof(RHS),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing})(::Array{Any,2}, ::Vararg{Any,N} where N) at /Users/dc/.julia/packages/DiffEqBase/1dZDc/src/diffeqfunction.jl:248
 [4] modelingtoolkitize(::ODEProblem{Array{Int64,2},Tuple{Float64,Float64},true,DiffEqBase.NullParameters,ODEFunction{true,typeof(RHS),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}) at /Users/dc/.julia/packages/ModelingToolkit/wy693/src/systems/diffeqs/modelingtoolkitize.jl:19
 [5] __solve(::ODEProblem{Array{Int64,2},Tuple{Float64,Float64},true,DiffEqBase.NullParameters,ODEFunction{true,typeof(RHS),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::MATLABDiffEq.ode15s, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}; saveat::Float64, timeseries_errors::Bool, reltol::Float64, abstol::Float64, kwargs::Base.Iterators.Pairs{Symbol,VectorContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing,Int64},Tuple{Symbol},NamedTuple{(:callback,),Tuple{VectorContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing,Int64}}}}) at /Users/dc/.julia/packages/MATLABDiffEq/b8nHE/src/MATLABDiffEq.jl:47
 [6] solve_call(::ODEProblem{Array{Int64,2},Tuple{Float64,Float64},true,DiffEqBase.NullParameters,ODEFunction{true,typeof(RHS),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::MATLABDiffEq.ode15s; merge_callbacks::Bool, kwargs::Base.Iterators.Pairs{Symbol,Any,Tuple{Symbol,Symbol},NamedTuple{(:callback, :saveat),Tuple{VectorContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing,Int64},Float64}}}) at /Users/dc/.julia/packages/DiffEqBase/1dZDc/src/solve.jl:96
 [7] #solve_up#454 at /Users/dc/.julia/packages/DiffEqBase/1dZDc/src/solve.jl:122 [inlined]
 [8] #solve#453 at /Users/dc/.julia/packages/DiffEqBase/1dZDc/src/solve.jl:106 [inlined]
 [9] top-level scope at /Users/dc/desktop/is/post.jl:106
 [10] include(::Module, ::String) at ./Base.jl:377
 [11] exec_options(::Base.JLOptions) at ./client.jl:288
 [12] _start() at ./client.jl:484
in expression starting at /Users/dc/desktop/is/post.jl:106

Thanks for any help!

2

There are a number of issues with your model. First off, MATLABDiffEq doesn’t allow for callbacks: it’ll throw an error when you get there anyways, so just use the normal methods (we cannot inject enough code into MATLAB to do it).

But secondly, the way you were referencing parameters wasn’t quite right. I think it might help if you use a LabelledArray, since then you can reference things by name. So I got a working code for you here:

using DifferentialEquations, LinearAlgebra, LabelledArrays

par = LVector(
    a = 25,
    d1 = 5,
    T = 2323,
    V = 1,
    c = 0.1,
    d2 = 0.35,
    e = 0.5,
    f = 0.2,
    g = 1.0/3.0,
    h = 0.3,
    k = 0.12,
    l = 1.0/11,
    m = 0.6,
    n = 0.1,
    p = n*l,
    q = (1-n)*l,
    z1 = 0.01,
    s = 0.2,
    t = q,
    u = 0.5*l,
    v = 0.592,
    w = v*u,
    x = (1-v)*u,
    y = 0.23245,
    z2  = 13,
    v1 = 0.0,
    v2 = 0.0)

function RHS(ydot,y,par,t)

    a = par.a;
    d = par[2];
    T = par.T;
    d = par[4];
    c = par.c;
    e = par.e;
    f = par.f;
    g = par.g;
    h = par.h;
    k = par.k;
    l = par.l;
    m = par.m;
    n = par.n;
    p = par.p;
    q = par.q;
    z = par[16];
    s = par.s;
    t = par.t;
    u = par.u;
    w = par.w;
    x = par.x;
    #y = par[22];
    z = par[23];
    v1= par.v1;
    v2= par.v2;

    fM = d * y[1] * (y[19] + e * (y[5] + y[8] + y[11] + y[14]) + f * y[17])/T;
    h = c * d * y[2] * (y[19] + e * (y[5] + y[8] + y[11] + y[14]) + f * y[17])/T;

    ydot[1] = -fM - (v1*y[1]-v2*y[2]);
    ydot[2] = -h + (v1*y[1]-v2*y[2]);
    ydot[3] = fM - k*y[3]-(v1*y[3]-v2*y[4]);
    ydot[4] = h - k*y[4]+(v1*y[3]-v2*y[4]);
    ydot[5] = k*(y[3]-y[5])-z*y[5]-(v1*y[5]-v2*y[6]);
    ydot[6] = k*(y[4]-y[6])-z*y[6]+(v1*y[5]-v2*y[6]);
    ydot[7] = z*(y[5]+y[6])-k*y[7];
    ydot[8] = k*(y[5]-y[8])-z*y[8]-(v1*y[8]-v2*y[9]);
    ydot[9] = k*(y[6]-y[9])-z*y[9]+(v1*y[8]-v2*y[9]);
    ydot[10] = z*(y[8]+y[9])-k*y[10] + k*y[7];
    ydot[11] = k*(y[8]-y[11])-z*y[11]-(v1*y[11]-v2*y[12]);
    ydot[12] = k*(y[9]-y[12])-z*y[12]+(v1*y[11]-v2*y[12]);
    ydot[13] = z*(y[11]+y[12])-k*y[13] + k*y[10];
    ydot[14] = k*(y[11]-y[14])-z*y[14]-(v1*y[14]-v2*y[15]);
    ydot[15] = k*(y[12]-y[15])-z*y[15]+(v1*y[14]-v2*y[15]);
    ydot[16] = z*(y[14]+y[15])-k*y[16] + k*y[13];
    ydot[17] = m*k*y[14]-z*y[17]-(v1*y[17]-v2*y[18])-t*y[17];
    ydot[18] = m*k*y[15]-z*y[18]+(v1*y[17]-v2*y[18])-t*y[18];
    ydot[19] = (1-m)*k*y[14]-(z+s)*y[19]-(v1*y[19]-v2*y[20])-(p+q)*y[19];
    ydot[20] = (1-m)*k*y[15]-(z+s)*y[20]+(v1*y[19]-v2*y[20])-(p+q)*y[20];
    ydot[21] = (1-m)*k*(y[14]+y[15]);

    return ydot
end

function condition(out,y,t,integrator) # Event when event_f(u,t) == 0
    out[1] = integrator.t - integrator.p[1] - integrator.p[2];
    out[2] = integrator.p[1] - integrator.t;
end

function affect!(integrator, idx)
    if idx == 1
        integrator.p.v1 = integrator.p[8]*integrator.p[9];
        integrator.p.v2 = integrator.p[8]*(1-integrator.p[9]);
    else
        integrator.p .= 0
    end
end

cb = VectorContinuousCallback(condition,affect!,2);
y0 = [2323.0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1];
tspan = (0.0, 100.0);
prob = ODEProblem(RHS,y0,tspan,par);
sol = solve(prob, Rodas5(), callback = cb, saveat = 1.0)

However, note that there are funny things in your model. You define d and z twice, so in the LabelledArray it caught that and threw an error. To get around that, I just renamed them d1 and d2, and z1 and z2, and left the way you were dereferencing in the ODE the same so you can see the oddities there (you’re grabbing it from an index for d1, so the other definition isn’t used? and why are some parameters like y non-existent?). So you probably need to clean up the definition of the ODE before it’s anything sensible, but in terms of the part for the question, this shows a clean way how to handle the parameter changes since then it’s just integrator.p.v1 = ... or to zero all parameters it’s just integrator.p .= 0.

1 Like
3

Thank you - that was very helpful! After I cleaned up my code, it runs properly now.

Since this was a minimal working example, I want to run my (highly) stiff system of ODEs across a wide range of parameter values. I say “highly” because of the discontinuities accounted for by the Callback functions.

As I’m running the solver in a loop over the parameter values, I noticed that the code freezes at a certain iteration (say, iteration 493 out of 100,000). This is after I supplied the full Jacobian and used Rodas5(). When I use Rodas4(), the code freezes at iteration 350.

So, my question is: I’m not sure which solver I should use, and what the values of reltol and abstol should be. Is there a general procedure to go about this?

Once again, thank you for your help!

4

Guess and check is usually pretty good. What do you mean by freezes? If you put an @show t in your function, is it getting seemingly stuck at a point? Is there something special at that point, like a discontinuity? Decreasing tolerances, like abstol=1e-8,reltol=1e-8 can increase the stability of solving some equations so it’s worth a try.

5

I see - yes, the solver is getting stuck at a point, so I’ll look into it.

6

In debugging, I noticed the following:

  1. “ydot” has some negative values at various time steps.
  2. The ODE integrator is stepping through the full time interval that was specified, for “proper” sets of parameter values.
  3. I might have specified the callbacks incorrectly…
  4. I would like to view the determinant and spectral radius of the Jacobian at each iteration. How can I do this?

For #3, my code is below:

function condition(out,y,t,integrator)
    out[1] = integrator.t - integrator.p[1] - integrator.p[2];
    out[2] = integrator.p[1] - integrator.t;
end

function affect!(integrator, idx)
    if idx == 1
        integrator.p[24] = integrator.p[8]*integrator.p[9];
        integrator.p[25] = integrator.p[8]*(1.0-integrator.p[9]);
    else
        integrator.y .= 0
        integrator.jacobian .= 0
    end
end

I believe the issue is in the else clause…

7

There are a few things in your equations that probably need fixing that I pointed out. Same variables are overwritten to multiple values, some aren’t used, etc. and I assume those are all accidents.

8

That’s not a thing. I don’t know what documentation you’re reading to come up with that, and I don’t know what you’re trying to do there.

9

Got it. Thanks for the feedback.

I’ve since figured out that the issue in my overall code is actually… not with the numerical integration routine. It’s elsewhere, and more manageable. It should be working soon.