Error solving ODE - Integrator stepped past tstops

Specific Domains Modelling & Simulations
1

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 : ERROR: LoadError: Something went wrong. Integrator stepped past tstops but the algorithm was dtchangeable. Please report this error. · GitHub
I’ve tried to kept it as simple as i could.

2

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.

3

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.

4

Fixed in OrdinaryDiffEq 3.21.0

3 Likes
5

Thank you for the quick answer and fix.