SciML/DiffEqFlux.jl

Error while training NeuralODE with exogenous inputs

mariaade26 opened this issue · 6 comments

mariaade26 commented

Hello everyone,

I'm currently in the process of training a NeuralODE model using data obtained from an ordinary differential equation (ODE) that describes the underlying physical phenomenon. In my code, I'm incorporating external disturbance data, Text and phih.

I've encountered an issue while working with two distinct types of disturbances represented by the variable phi. One of these disturbances seems to be causing problems, while the other one works as expected. In one version of my code, I'm able to successfully train the neural network, but in the other, I encounter issues such as DimensionMismatch errors and instability warnings during the simulation.
I'm seeking assistance in understanding the root cause of this problem. Can anyone provide insights or suggestions?

WORKING CODE

using DifferentialEquations, Lux, ComponentArrays, DiffEqFlux, Optimization,  OptimizationOptimJL, Plots, Random, OptimizationOptimisers, DataInterpolations, CSV

text = open(f->read(f, String), "File csv/eplusout.csv") |> s -> split(s, "\n")[1:end-1]
data=parse.(Float32, text)  #valori di temperatura esterna ogni minuto per una settimana circa
Testerna= data[1:1800];

phi1=0*ones(300)
phi2=5*ones(300)
phi3=10*ones(300)
phi4=30*ones(300)
phi5=20*ones(300)
phi6=0*ones(300)

phih=vcat(Float32.(phi1),Float32.(phi2),Float32.(phi3),Float32.(phi4),Float32.(phi5),Float32.(phi6));
tspan= (0.0f0, 1800.0f00)
t= range(tspan[1], tspan[2], length=1800)

disturbance = LinearInterpolation(Testerna,t);
disturbance2 = LinearInterpolation(phih,t);

function Text(t)
    return disturbance(t)
end

function phi(t)
    return disturbance2(t)
end 

function f(du, u, p, t)
    Rv = p[1]
    Ci = p[2]
    du[1] = 1/(Rv*Ci) .* (Text(t) .- u[1]) .+ phi(t)/Ci

end 
 
u0= [5.0]
p= [0.05, 100.0]


prob = ODEProblem(f, u0, tspan, p)
sol = solve(prob, Tsit5(), saveat=t)
tplot=t
plot(sol, xlabel="Tempo", ylabel="Temperatura [°C]", legend=false)

rng=Random.default_rng()

nn_model = Lux.Chain(Lux.Dense(3, 8, tanh), Lux.Dense(8, 1))
p_model, st = Lux.setup(rng, nn_model)

function dudt(u, p, t)
    global st
    out, st = nn_model(vcat(u[1], Text(t),phi(t)), p, st)
    return out
end

prob = ODEProblem(dudt, u0, tspan, nothing)

function predict_neuralode(p)
    _prob = remake(prob, p = p)
    Array(solve(_prob, Tsit5(), saveat = t, abstol = 1e-8, reltol = 1e-6))
end

function loss(p)
    pred = predict_neuralode(p)
    N = length(pred)
    return sum(abs2.(sol[1,:] .- pred')) 
end

callback = function(p,l)
        println(l)
     return false 
    end 

adtype = Optimization.AutoZygote()
pinit= ComponentArray(p_model)
optf = Optimization.OptimizationFunction((x, p) -> loss(x), adtype)
optprob = Optimization.OptimizationProblem(optf, pinit)
@time begin 
res0 = Optimization.solve(optprob, Optimisers.Adam(0.05), callback=callback, maxiters = 300)
end 

pred= predict_neuralode(res0.u)
pred=pred'
plot(t,pred, label= "prediction")
plot!(t, sol[1,:], label= "data", xlabel="Tempo", ylabel="Temperatura [°C]") 

optprob2 = remake(optprob,u0 = res0.u)

@time begin
result_neuralode2 = Optimization.solve(optprob2,
                                        Optim.BFGS(initial_stepnorm=0.01),
                                        callback=callback,
                                        allow_f_increases = false)
end 
pred2= predict_neuralode(result_neuralode2.u)
pred2=pred2'
plot(t,pred2, label= "prediction")
plot!(t, sol[1,:], label= "data", xlabel="Tempo", ylabel="Temperatura [°C]")

In the non-working version, I only replaced these lines

"phi1=0ones(300)
phi2=5ones(300)
phi3=10ones(300)
phi4=30ones(300)
phi5=20ones(300)
phi6=0ones(300)
phih=vcat(Float32.(phi1),Float32.(phi2),Float32.(phi3),Float32.(phi4),Float32.(phi5),Float32.(phi6));"

with the following


 text = open(f->read(f, String), "File csv/phi heating .csv") |> s -> split(s, "\n")[1:end-1]
phih=parse.(Float32, text);

In particular, as soon as I run the simulation with this modification I get this errore "┌ Warning: dt(0.00012207031) <= dtmin(0.00012207031) at t=179.09944, and step error estimate = 1.4775506934013858. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase C:\Users\Michele.julia\packages\SciMLBase\kTUaf\src\integrator_interface.jl:599
DimensionMismatch: arrays could not be broadcast to a common size; got a dimension with lengths 1800 and 179".

mariaade26 commented

NON WORKING CODE

using DifferentialEquations, Lux, ComponentArrays, DiffEqFlux, Optimization,  OptimizationOptimJL, Plots, Random, OptimizationOptimisers, DataInterpolations, CSV

text = open(f->read(f, String), "File csv/eplusout.csv") |> s -> split(s, "\n")[1:end-1]
data=parse.(Float32, text)  #valori di temperatura esterna ogni minuto per una settimana circa
Testerna= data[1:1800];

text = openphih=parse.(Float32, text);(f->read(f, String), "File csv/phi heating .csv") |> s -> split(s, "\n")[1:end-1]
 phih=parse.(Float32, text);

tspan= (0.0f0, 1800.0f00)
t= range(tspan[1], tspan[2], length=1800)

disturbance = LinearInterpolation(Testerna,t);
disturbance2 = LinearInterpolation(phih,t);

function Text(t)
    return disturbance(t)
end

function phi(t)
    return disturbance2(t)
end 

function f(du, u, p, t)
    Rv = p[1]
    Ci = p[2]
    du[1] = 1/(Rv*Ci) .* (Text(t) .- u[1]) .+ phi(t)/Ci

end 
 
u0= [5.0]
p= [0.05, 100.0]


prob = ODEProblem(f, u0, tspan, p)
sol = solve(prob, Tsit5(), saveat=t)
tplot=t
plot(sol, xlabel="Tempo", ylabel="Temperatura [°C]", legend=false)

rng=Random.default_rng()

nn_model = Lux.Chain(Lux.Dense(3, 8, tanh), Lux.Dense(8, 1))
p_model, st = Lux.setup(rng, nn_model)

function dudt(u, p, t)
    global st
    out, st = nn_model(vcat(u[1], Text(t),phi(t)), p, st)
    return out
end

prob = ODEProblem(dudt, u0, tspan, nothing)

function predict_neuralode(p)
    _prob = remake(prob, p = p)
    Array(solve(_prob, Tsit5(), saveat = t, abstol = 1e-8, reltol = 1e-6))
end

function loss(p)
    pred = predict_neuralode(p)
    N = length(pred)
    return sum(abs2.(sol[1,:] .- pred')) 
end

callback = function(p,l)
        println(l)
     return false 
    end 

adtype = Optimization.AutoZygote()
pinit= ComponentArray(p_model)
optf = Optimization.OptimizationFunction((x, p) -> loss(x), adtype)
optprob = Optimization.OptimizationProblem(optf, pinit)
@time begin 
res0 = Optimization.solve(optprob, Optimisers.Adam(0.05), callback=callback, maxiters = 300)
end 

pred= predict_neuralode(res0.u)
pred=pred'
plot(t,pred, label= "prediction")
plot!(t, sol[1,:], label= "data", xlabel="Tempo", ylabel="Temperatura [°C]") 

optprob2 = remake(optprob,u0 = res0.u)

@time begin
result_neuralode2 = Optimization.solve(optprob2,
                                        Optim.BFGS(initial_stepnorm=0.01),
                                        callback=callback,
                                        allow_f_increases = false)
end 
pred2= predict_neuralode(result_neuralode2.u)
pred2=pred2'
plot(t,pred2, label= "prediction")
plot!(t, sol[1,:], label= "data", xlabel="Tempo", ylabel="Temperatura [°C]")
prbzrg commented

In particular, as soon as I run the simulation with this modification I get this errore "┌ Warning: dt(0.00012207031) <= dtmin(0.00012207031) at t=179.09944, and step error estimate = 1.4775506934013858. Aborting. There is either an error in your model specification or the true solution is unstable.

It might help to remove Tsit5() from those two solve calls.

DimensionMismatch: arrays could not be broadcast to a common size; got a dimension with lengths 1800 and 179".

It might help to read data like:

using CSV, DataFrames, MLJModelInterface

fn = raw"csv/phi heating.csv"
phih_df = CSV.read(fn, DataFrame) # 👈 read as DataFrame

# any manipulation of data 👇

#######################

phih = collect(transpose(MLJModelInterface.matrix(phih_df)))  # 👈 convert it to an array (each column is a data point)
mariaade26 commented

Thank you a lot for your help, but that didn't solve my problem!
I read the csv as you said and instead of Tsit5() I used AutoTsit5(Rosenbrock23())

mariaade26 commented

In particular, as soon as I run the simulation with this modification I get this errore "┌ Warning: dt(0.00012207031) <= dtmin(0.00012207031) at t=179.09944, and step error estimate = 1.4775506934013858. Aborting. There is either an error in your model specification or the true solution is unstable.

It might help to remove Tsit5() from those two solve calls.

DimensionMismatch: arrays could not be broadcast to a common size; got a dimension with lengths 1800 and 179".

It might help to read data like:

using CSV, DataFrames, MLJModelInterface

fn = raw"csv/phi heating.csv"
phih_df = CSV.read(fn, DataFrame) # 👈 read as DataFrame

# any manipulation of data 👇

#######################

phih = collect(transpose(MLJModelInterface.matrix(phih_df)))  # 👈 convert it to an array (each column is a data point)

@prbzrg if you want to try to run it, ifever you had time, here is the csv files.
phi heating .csv
Temperatura esterna.csv

prbzrg commented

Changing the algorithm in solve and reltol and abstol values fixed the problem for me.
Try this:

using DifferentialEquations,
    Lux,
    ComponentArrays,
    DiffEqFlux,
    Optimization,
    OptimizationOptimJL,
    Plots,
    Random,
    OptimizationOptimisers,
    DataInterpolations,
    CSV,
    DataFrames

data_fn = raw"csv/Temperatura.esterna.csv"
data_df = CSV.read(data_fn, DataFrame; header = false);

Testerna = convert.(Float32, data_df[1:1800, 1]);

phih_fn = raw"csv/phi.heating.csv"
phih_df = CSV.read(phih_fn, DataFrame; header = false);

phih = convert.(Float32, phih_df[!, 1]);

tspan = (0.0f0, 1800.0f0)
t = range(tspan[1], tspan[2]; length = 1800)

disturbance = LinearInterpolation(Testerna, t);
disturbance2 = LinearInterpolation(phih, t);

function Text(t)
    return disturbance(t)
end

function phi(t)
    return disturbance2(t)
end

function f(du, u, p, t)
    Rv = p[1]
    Ci = p[2]
    du[1] = 1.0f0 / (Rv * Ci) .* (Text(t) .- u[1]) .+ phi(t) / Ci
end

u0 = [5.0f0]
p = [0.05f0, 100.0f0]

prob = ODEProblem(f, u0, tspan, p)
sol = solve(
    prob,
    VCABM();
    saveat = t,
    reltol = sqrt(eps(one(Float32))),
    abstol = eps(one(Float32)),
)  # 👈 changed
tplot = t
plot(sol; xlabel = "Tempo", ylabel = "Temperatura [°C]", legend = false)

rng = Random.default_rng()

nn_model = Lux.Chain(Lux.Dense(3, 8, tanh), Lux.Dense(8, 1))
p_model, st = Lux.setup(rng, nn_model)

function dudt(u, p, t)
    global st
    out, st = nn_model(vcat(u[1], Text(t), phi(t)), p, st)
    return out
end

prob = ODEProblem(dudt, u0, tspan, nothing)

function predict_neuralode(p)
    _prob = remake(prob; p = p)
    Array(
        solve(
            _prob,
            VCABM();
            saveat = t,
            reltol = sqrt(eps(one(Float32))),
            abstol = eps(one(Float32)),
        ),
    )  # 👈 changed
end

function loss(p)
    pred = predict_neuralode(p)
    N = length(pred)
    return sum(abs2.(sol[1, :] .- pred'))
end

callback = function (p, l)
    println(l)
    return false
end

adtype = Optimization.AutoZygote()
pinit = ComponentArray(p_model)
optf = Optimization.OptimizationFunction((x, p) -> loss(x), adtype)
optprob = Optimization.OptimizationProblem(optf, pinit)
@time begin
    res0 = Optimization.solve(
        optprob,
        Optimisers.Adam(0.05f0);
        callback = callback,
        maxiters = 300,
    )
end

pred = predict_neuralode(res0.u)
pred = pred'
plot(t, pred; label = "prediction")
plot!(t, sol[1, :]; label = "data", xlabel = "Tempo", ylabel = "Temperatura [°C]")

optprob2 = remake(optprob; u0 = res0.u)

@time begin
    result_neuralode2 = Optimization.solve(
        optprob2,
        Optim.BFGS(; initial_stepnorm = 0.01f0);
        callback = callback,
        allow_f_increases = false,
    )
end
pred2 = predict_neuralode(result_neuralode2.u)
pred2 = pred2'
plot(t, pred2; label = "prediction")
plot!(t, sol[1, :]; label = "data", xlabel = "Tempo", ylabel = "Temperatura [°C]")
mariaade26 commented

Thank you very very much!
I've just tried this and it seems to work!
I'm immensely grateful