Error solving ODE - Integrator stepped past tstops

I get a error when solving a differential equation with the DifferentialEquations package. The error appear or not depending on the seed i’m specifying. I can’t understand from where this error come from, i’m not sure if it’s related to something i’m doing wrong or related to the package, and where i should post this issue.


ERROR: LoadError: Something went wrong. Integrator stepped past tstops but the algorithm was dtchangeable. Please report this error.

Stacktrace:
 [1] handle_tstop!(::OrdinaryDiffEq.ODEIntegrator{OrdinaryDiffEq.Tsit5,Array{Float64,2},Float64,Tuple{Array{Float64,2},Array{Float64,2},Array{Float64,2},Int64,Int64,Array{Any,1},Array{Any,1},Array{Float64,1},Float64,Float64},Float64,Float64,Float64,Array{Array{Float64,2},1},DiffEqBase.ODESolution{Float64,3,Array{Array{Float64,2},1},Void,Void,Array{Float64,1},Array{Array{Array{Float64,2},1},1},DiffEqBase.ODEProblem{Array{Float64,2},Float64,true,Tuple{Array{Float64,2},Array{Float64,2},Array{Float64,2},Int64,Int64,Array{Any,1},Array{Any,1},Array{Float64,1},Float64,Float64},#water_heater_dyn_classique_2layers!,Void,Void,UniformScaling{Int64},DiffEqBase.StandardODEProblem},OrdinaryDiffEq.Tsit5,OrdinaryDiffEq.InterpolationData{#water_heater_dyn_classique_2layers!,Array{Array{Float64,2},1},Array{Float64,1},Array{Array{Array{Float64,2},1},1},OrdinaryDiffEq.Tsit5Cache{Array{Float64,2},Array{Float64,2},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}}}},#water_heater_dyn_classique_2layers!,Void,OrdinaryDiffEq.Tsit5Cache{Array{Float64,2},Array{Float64,2},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,DiffEqBase.#ODE_DEFAULT_NORM,DiffEqBase.CallbackSet{Tuple{},Tuple{DiffEqBase.DiscreteCallback{#condition2#59,#affect2!#60{Array{Int64,1}},DiffEqBase.#INITIALIZE_DEFAULT},DiffEqBase.DiscreteCallback{#condition3#61,#affect3!#62,DiffEqBase.#INITIALIZE_DEFAULT},DiffEqBase.DiscreteCallback{#condition4#63{Float64},#affect4!#64,DiffEqBase.#INITIALIZE_DEFAULT},DiffEqBase.DiscreteCallback{#condition5#65{Float64},#affect5!#66,DiffEqBase.#INITIALIZE_DEFAULT},DiffEqBase.DiscreteCallback{#condition6#67{Float64},#affect6!#68,DiffEqBase.#INITIALIZE_DEFAULT},DiffEqBase.DiscreteCallback{#condition7#69{Float64},#affect7!#70,DiffEqBase.#INITIALIZE_DEFAULT},DiffEqBase.DiscreteCallback{DiffEqJump.DirectJumpAggregation{Float64,DiffEqJump.MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{#rate#57{Array{Float64,2},Int64}},Tuple{#affect!#58{Int64}},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,DiffEqJump.MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{#rate#57{Array{Float64,2},Int64}},Tuple{#affect!#58{Int64}},RandomNumbers.Xorshifts.Xoroshiro128Star},DiffEqJump.DirectJumpAggregation{Float64,DiffEqJump.MassActionJump{Array{Float64,1},Array{Array{Pair{Int64,Float64},1},1},Array{Array{Pair{Int64,Float64},1},1}},Tuple{#rate#57{Array{Float64,2},Int64}},Tuple{#affect!#58{Int64}},RandomNumbers.Xorshifts.Xoroshiro128Star}}}},DiffEqBase.#ODE_DEFAULT_ISOUTOFDOMAIN,DiffEqBase.#ODE_DEFAULT_PROG_MESSAGE,DiffEqBase.#ODE_DEFAULT_UNSTABLE_CHECK,DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Void,Void,Int64,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Array{Float64,2}}) at /home/piedmari/.julia/v0.6/OrdinaryDiffEq/src/integrators/integrator_utils.jl:417

Here is a example of a code where it happen : https://gist.github.com/mariepied/2547a1a740bb76e57d979a29cde7ffd1
I’ve tried to kept it as simple as i could.

That’s quite an example. I won’t be able to get to it until maybe tonight or tomorrow on a plane. Since I won’t have my whole setup, it would be helpful if we can narrow this down a little bit more. Does it need for k in 1:Nb_wh or is there one solve in there that is an issue? Does it need every callback to hit the problem? If you made the ODE part trivial (dx .= 0) do you still hit the issue?

There’s a chance it’s a floating point error that arises only when you put everything together, but that should be unlikely.

Ugh it’s due to the fact that floating point isn’t associative, so its off by a little bit:

(top(tstops), integrator.t) = (743.8332651649943, 117.37162136488763)
top(tstops) - (integrator.t + nextfloat(integrator.dt)) = -1.1368683772161603e-13
integrator.dt = 626.4616438001067
integrator.t = 743.8332651649944

I have a fix in mind. More minimal example:

import DifferentialEquations
import Interpolations
import Distributions
using Plots

function module_simulation_wh(CONT,trajectory,mu_opt,wh,x_init,theta_0,heat_on)
    Nb_wh,n_l,m_c,control_pos,K,maxPower,xenv,xL,densW,V,M,Cpf,U,Ai,R,qx0,L,delta_t,T,ft,n_t,y,lowComfort,highComfort,B,A0,A1,H,zeta0 = wh

    x_tot = zeros(Nb_wh,n_l,n_t)
    theta_tot = zeros(Nb_wh,n_t)
    x_mean_tot = zeros(Nb_wh,n_t)
    u_tot = zeros(Nb_wh,m_c,n_t)
    uarray_tot = zeros(Nb_wh,n_t)

    index = 2

    for k in 1:Nb_wh
        @show k

        ## Markov Chain ##
        rate(u,p,t) = L[2,1]*u[1+index,1] +L[1,2]*(1-u[1+index,1]);

        function affect!(integrator)
            if integrator.u[1+index,1] == 0;
                integrator.u[1+index,:] = 1;
            else
                integrator.u[1+index,:] = 0;
            end
        end
        jump_0 = DifferentialEquations.JumpSet(DifferentialEquations.ConstantRateJump(rate,affect!));

        # Simulation #
        x = zeros(n_l,n_t);
        theta = zeros(n_t,);
        x[:,1] = x_init[k,:]

        x_mean = zeros(n_t,)
        theta_bis = repmat([theta_0[k]],1,n_l)
        x_0 = [x_init[k,:] heat_on[k,:] theta_bis']'
        u = zeros(m_c,n_t)
        uarray = zeros(n_t,)
        c0 = Array{Any}(1,)
        c1 = Array{Any}(1,)
        ufree = Array{Any}(1,)
        c0[1],c1[1],ufree[1] = calcul_cufree(n_l,M,B,U,Ai,mean(x_init[k,:]),xenv,K,Cpf,xL,my_zeta(T,L,zeta0))
        prob = DifferentialEquations.ODEProblem(water_heater_dyn_classique_2layers!,x_0,(0.0,T),(A0,A1,B,1,maxPower,c0,c1,H,highComfort,lowComfort))
        jump_prob = DifferentialEquations.JumpProblem(prob,DifferentialEquations.Direct(),jump_0);
        sol = DifferentialEquations.solve(jump_prob,DifferentialEquations.Tsit5());
    end
    return x_tot,u_tot,theta_tot,x_mean_tot,uarray_tot,trajectory,CONT

end


function my_zeta(t,L,zeta0)
    ## t time, L: infenitesimal generator (n_etat,n_etat), zeta0 (n_etat,)
    return zeta0'*expm(t*L)
end

function calcul_cufree(n,M,B,U,Ai,xinitmean,xenv,K,Cpf,xL,zeta_)
    temp = zeros(n,)
    temp[n] = 1
    ufree1,ufree2 = calcul_ufree(n,M,B,U,Ai,xinitmean,xenv,K,Cpf,xL)
    # Vector C / c0 if m_dot = 0, c1 si m_dot = K %
    c0=U*Ai/(M*Cpf)*xenv
    c1=c0+K/M*xL*temp
    return c0,c1,ufree1+zeta_[2]*ufree2
end

function calcul_ufree(n,M,B,U,Ai,xinitmean,xenv,K,Cpf,xL)
    ufree_1 = sum(U*Ai.*(xinitmean-xenv))/n*ones(n,)/hour_seconds
    ufree_2 = zeros(n,)
    ufree_2[n] = K*Cpf*(xinitmean[1]-xL)/hour_seconds
    return ufree_1,ufree_2
end


function water_heater_dyn_classique_2layers!(dx,x,p,t)
    ## x = (temperature , heating state, etat(0,1) of the markov chain)
    A0,A1,B,N_wh,maxPower,c0,c1,H,highComfort,lowComfort = p
    # c0,c1 : Array{Array(n,)}(N_wh,1)
    # sol_pi : Array{solution Diffeq}(N_wh,1)
    # sol_offset : Array{solution Diffeq}(N_wh,1)
    # ufree : Array{Array(m,)}(N_wh,1)

    m = size(B)[2]
    len = length(H)
    for k in 1:N_wh
        u = x[k+N_wh,:]*maxPower/2
        if x[k+2*N_wh,1] == 0
            dx[k,:] = A0*x[k,:]+B*u+c0[k];
        else
            dx[k,:] = A1*x[k,:]+B*u+c1[k];
        end
        dx[k+N_wh,:] = 0
        dx[k+2*N_wh,:] = 0
    end
end

hour_seconds = 3600
minute_seconds = 60
Kelvin = 273.15


Nb_wh=1100
n_l=2
m_c=2
control_pos= [1,n_l]
K=1.31*60/hour_seconds
maxPower=4500
xenv=25+Kelvin
xL= 15+Kelvin
densW=1000
V=0.2730
M = V*densW/n_l
Cpf=4190
U=28.38*60/hour_seconds
Ai = [1.27825, 1.27825]
R=0.025*eye(n_l)/hour_seconds
qx0=8000/hour_seconds
L=[[-0.5 0.5];[7 -7]]/hour_seconds
lowComfort=50+Kelvin
highComfort=60+Kelvin
B = [1.74845e-6 0.0; 0.0 1.74845e-6]
A0 = [-1.05713e-6 -0.0; -0.0 -1.05713e-6]
A1 = [-0.000161008 0.000159951; 0.0 -0.000161008]
H= 1/n_l*ones(n_l,)

zeta0 = [0.5; 0.5]

### Simulation ####
srand(1819275)
delta_t = minute_seconds
T= 2.0*hour_seconds
ft = 0:delta_t:T
n_t = length(ft)
xinitmean = (53+Kelvin)*ones(2,)
y=55+Kelvin
z = lowComfort

wh = (Nb_wh,n_l,m_c,control_pos,K,maxPower,xenv,xL,densW,V,M,Cpf,U,Ai,R,qx0,L,delta_t,T,ft,n_t,y,lowComfort,highComfort,B,A0,A1,H,zeta0)

x_init = zeros(Nb_wh,n_l)
for j in 1:n_l
    x_init[:,j] = rand(Distributions.Normal(xinitmean[j],1),Nb_wh);
    for k in 1:Nb_wh
        if j>=2 && x_init[k,j] > x_init[k,j-1]
            temp = x_init[k,j-1]
            x_init[k,j-1] = x_init[k,j]
            x_init[k,j] = temp
        end
    end
end
theta_0 = rand([0 1],Nb_wh)
heat_on = ones(Nb_wh,)
heat_on[1:round(Int,Nb_wh/3*2)] = 0
heat_on = shuffle(heat_on)
heat_on = repmat(heat_on,1,n_l)
heat_on[:,1] = zeros(Nb_wh)

x,u,theta,x_mean,uarray,trajectory2,CONT = module_simulation_wh(false,[],0,wh,x_init,theta_0,heat_on)

It’s hard to get more precise than that because it’s due to the scaling of eps in floating point calculations. I’ll try a fix in a sec.

Fixed in OrdinaryDiffEq 3.21.0

3 Likes

Thank you for the quick answer and fix.