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

error_diffEq.jl

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 
        
        ## 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!));

            # Callbacks #
            cb = DifferentialEquations.CallbackSet()
            for i in 2:m_c
                for j in control_pos[i]:-1:2
                    
                    function condition2(u,t,integrator)
                        (u[2,i] == 1 && u[1,j] > u[1,j-1])
                    end

                    function affect2!(integrator)
                        integrator.u[1,j-1:control_pos[i]] = mean(integrator.u[1,j-1:control_pos[i]])
                    end
                    cb = DifferentialEquations.CallbackSet(cb,DifferentialEquations.DiscreteCallback(condition2,affect2!,save_positions=(true,true)))

                end
            end
    
            for l in n_l:-1:2
                for p in l:-1:2
                    function condition3(u,t,integrator)
                        (u[1,l] > u[1,p-1])
                    end

                    function affect3!(integrator)
                        integrator.u[1,p-1:l] = mean(integrator.u[1,p-1:l])
                    end
                    cb = DifferentialEquations.CallbackSet(cb,DifferentialEquations.DiscreteCallback(condition3,affect3!,save_positions=(true,true)))

                end
            end


            function condition4(u,t,integrator)
                (u[2,1] == 1 && u[1,1] >= highComfort)
            end

            function affect4!(integrator)
                integrator.u[2,1] = 0
            end
            cb = DifferentialEquations.CallbackSet(cb,DifferentialEquations.DiscreteCallback(condition4,affect4!,save_positions=(true,true)))

            function condition5(u,t,integrator)
                (u[2,1] == 0 && u[1,1] < lowComfort)
            end

            function affect5!(integrator)
                integrator.u[2,1] = 1
                integrator.u[2,2] = 0
            end
            cb = DifferentialEquations.CallbackSet(cb,DifferentialEquations.DiscreteCallback(condition5,affect5!,save_positions=(true,true)))

            function condition6(u,t,integrator)
                (u[2,2] == 1 && u[1,2] >= highComfort)
            end

            function affect6!(integrator)
                integrator.u[2,2] = 0
            end
            cb = DifferentialEquations.CallbackSet(cb,DifferentialEquations.DiscreteCallback(condition6,affect6!,save_positions=(true,true)))

            function condition7(u,t,integrator)
                (u[2,2] == 0 && u[1,2] < lowComfort && u[2,1] == 0)
            end

            function affect7!(integrator)
                integrator.u[2,2] = 1
            end
            cb = DifferentialEquations.CallbackSet(cb,DifferentialEquations.DiscreteCallback(condition7,affect7!,save_positions=(true,true)))
            
            # 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,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(),callback = cb,dtmax=0.01*hour_seconds);
            len = length(H)
            for i in 1:n_l
                x[i,:] = [sol(temp)[1,i] for temp in ft];
            end
            for i in 1:n_t
                u[:,i] = sol(ft[i])[2,:]
            end
            theta = [sol(temp)[3,1] for temp in ft];
            x_mean = mean(x,1)
            uarray = sum(u,1)

        x_tot[k,:,:] = x 
        theta_tot[k,:] = theta
        x_mean_tot[k,:] = x_mean
        u_tot[k,:,:] = u
        uarray_tot[k,:] = uarray
    end
    return x_tot,u_tot,theta_tot,x_mean_tot,uarray_tot,trajectory,CONT

end


function 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)


plotly()
plot(ft/hour_seconds,x_mean[1:10,:]'-Kelvin)
display(plot!(ft/hour_seconds,mean(x_mean[:,:],1)'-Kelvin,linecolor=:red,title = "Temperature evolution in 10 water-heaters.",xlabel = "Temps en h",ylabel = "°C"))
display(plot(ft/hour_seconds,sum(uarray[:,:],1)'/Nb_wh,legend = false,title = "Water-Heater Energy consumption per house",xlabel = "Temps en h",ylabel = "watts"))