Training a NeuralODE with an ODE depending on exogenous time-dependent input

General Usage Performance
1

Hi everyone.
I’m working on developing a NeuralODE trained on ad ODE which describes a building’s thermal behaviour through an RC equation.

I tried both with the definition of a NeuralODE from a multilayer NN and by training a single layer NN with Flux.train!. In both cases I get a dimensional mismatch problem.
I followed some examples solved by @ChrisRackauckas such as How to create a multilayer perceptron using DiffEqFlux.jl · Issue #554 · SciML/DiffEqFlux.jl · GitHub, but still I couldn’t solve/understand this issues.

I would really appreciate some help!
Thanks in advace.

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

function RC!(du,u,p,t)
Rv, Ci = p
Text = ext_disturbance(t)
phih= ext_disturbance2(t)

du[1] = 1/(Rv*Ci) .* (Text .- u[1]) .+ phih/Ci

end

u0= [5.0]
tspan= (0.0f0, 5.0f00)
t= range(tspan[1], tspan[2], length=10)
p= [0.05, 10.0]

Testerna = [1,2,3,4,5,6,7,8,9,10]
phih = [10, 20, 30, 40, 50, 60 ,70 ,80 ,90 ,100];

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

function ext_disturbance(t)
return disturbance(t)
end

function ext_disturbance2(t)
return disturbance2(t)
end

ode

prob=ODEProblem(RC!, u0, tspan, p)
sol= solve(prob, Tsit5(), saveat=t)
plot(t,sol[1,:], label = “reale”,xlabel = “tempo [s]”, ylabel = “Temperatura [°C]”);

prob3 = ODEProblem(RC!, u0, tspan, p_random)
dudt= Lux.Chain( p-> solve(prob3, Tsit5(), saveat=t, p=p)[1,:],
vec,
Lux.Dense(10,50, tanh),
Lux.Dense(50,1))

n_ode= NeuralODE(dudt, tspan, Tsit5(), saveat=t)
p = [0.05, 10]
rng= Random.default_rng()
p, st = Lux.setup(rng,dudt)

function predict_neuralode(p)
Array(n_ode(u0,p,st)[1])
end

function loss_neuralode(p)
pred= predict_neuralode(p)
loss = sum(abs2, sol .- pred)
return loss, pred
end

callback = function ()
display(loss())
end

pinit = ComponentArray(p)

adtype = Optimization.AutoZygote()

optf = Optimization.OptimizationFunction((x, p) → loss_neuralode(x), adtype)
optprob = Optimization.OptimizationProblem(optf, pinit)

result_neuralode = Optimization.solve(optprob,
Optimisers.Adam(0.05),
callback = callback,
maxiters = 400)

This is the beginning of the error that I got ( it is too long to be copied):
BoundsError: attempt to access 1-element Vector{Float64} at index [2]

Stacktrace:
[1] getindex
@ .\essentials.jl:13 [inlined]
[2] indexed_iterate
@ .\tuple.jl:89 [inlined]
[3] RC!(du::Vector{Float64}, u::Vector{Float64}, p::Vector{Float64}, t::Float32)
@ Main .\In[78]:5
[4] Void
@ C:\Users\Michele.julia\packages\SciMLBase\s9wrq\src\utils.jl:473 [inlined]
[5] (::FunctionWrappers.CallWrapper{Nothing})(f::SciMLBase.Void{typeof(RC!)}, arg1::Vector{Float64}, arg2::Vector{Float64}, arg3::Vector{Float64}, arg4::Float32)
@ FunctionWrappers C:\Users\Michele.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:65
[6] macro expansion
@ C:\Users\Michele.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:137 [inlined]
[7] do_ccall
@ C:\Users\Michele.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:125 [inlined]
[8] FunctionWrapper
@ C:\Users\Michele.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:144 [inlined]
[9] _call
@ C:\Users\Michele.julia\packages\FunctionWrappersWrappers\9XR0m\src\FunctionWrappersWrappers.jl:12 [inlined]
[10] FunctionWrappersWrapper
@ C:\Users\Michele.julia\packages\FunctionWrappersWrappers\9XR0m\src\FunctionWrappersWrappers.jl:10 [inlined]
[11] ODEFunction
@ C:\Users\Michele.julia\packages\SciMLBase\s9wrq\src\scimlfunctions.jl:2126 [inlined] …

2

With Flux.train! the problem showed like this:

‘’’ p_random= [-2.0, 5.0]
prob2= ODEProblem(RC!, u0, tspan, p_random)
pred() = solve(prob2, Tsit5(), saveat=t)

flux TRAIN

param= Flux.params(p_random)
loss()=Flux.mse(pred()[1,:], sol[1,:])
data = Iterators.repeated((),100)
opt= ADAM(0.1)

cb = function()
display(loss())
plt = plot(t, sol, label = “data”)
plot!(plt, t, pred[1,:], label = “prediction”)

display(plot(plt))

end

Flux.train!(loss, param, data, opt, cb=cb) ‘’’

with the following error:
753.2408068488701
Expects 10 elements in each col of y, found 1.

Stacktrace:
[1] error(s::String)
@ Base .\error.jl:35
[2] _compute_xyz(x::StepRangeLen{Float32, Float64, Float64, Int64}, y::Vector{Float64}, z::Nothing, nice_error::Bool)
@ RecipesPipeline C:\Users\Michele.julia\packages\RecipesPipeline\BGM3l\src\series.jl:90
[3] macro expansion
@ C:\Users\Michele.julia\packages\RecipesPipeline\BGM3l\src\series.jl:140 [inlined]
[4] apply_recipe(plotattributes::AbstractDict{Symbol, Any}, #unused#::Type{RecipesPipeline.SliceIt}, x::Any, y::Any, z::Any)
@ RecipesPipeline C:\Users\Michele.julia\packages\RecipesBase\BRe07\src\RecipesBase.jl:300
[5] _process_userrecipes!(plt::Any, plotattributes::Any, args::Any)
@ RecipesPipeline C:\Users\Michele.julia\packages\RecipesPipeline\BGM3l\src\user_recipe.jl:38
[6] recipe_pipeline!(plt::Any, plotattributes::Any, args::Any)
@ RecipesPipeline C:\Users\Michele.julia\packages\RecipesPipeline\BGM3l\src\RecipesPipeline.jl:72
[7] _plot!(plt::Plots.Plot, plotattributes::Any, args::Any)
@ Plots C:\Users\Michele.julia\packages\Plots\rz1WP\src\plot.jl:223
[8] plot(::Any, ::Vararg{Any}; kw::Base.Pairs{Symbol, V, Tuple{Vararg{Symbol, N}}, NamedTuple{names, T}} where {V, N, names, T<:Tuple{Vararg{Any, N}}})
@ Plots C:\Users\Michele.julia\packages\Plots\rz1WP\src\plot.jl:102
[9] (::var"#60#61")()
@ Main .\In[85]:14
[10] macro expansion
@ C:\Users\Michele.julia\packages\Flux\n3cOc\src\optimise\train.jl:149 [inlined]
[11] macro expansion
@ C:\Users\Michele.julia\packages\ProgressLogging\6KXlp\src\ProgressLogging.jl:328 [inlined]
[12] train!(loss::Function, ps::Zygote.Params{Zygote.Buffer{Any, Vector{Any}}}, data::Base.Iterators.Take{Base.Iterators.Repeated{Tuple{}}}, opt::Adam; cb::var"#60#61")
@ Flux.Optimise C:\Users\Michele.julia\packages\Flux\n3cOc\src\optimise\train.jl:140
[13] top-level scope
@ In[85]:21

3

The answer is to use DataInterpolations to represent the time-dependent input via a spline. More in-depth response will come in the future but low on time right now.

4

Hi, thank you very much for your answer.
In the last days I managed to solve several issues that I had, but a problem still appears when I want to add the ODE solver into the neural network as a layer.
I couldn’t fix it, even by looking at the available examples, therefore I would be glad if you could help me!
Here is the code that I managed to use by adding (x → x) to the NN, instead of my solver, but this is not good when the ODE gets more complicated.
Thank you in advance.


function RC!(du,u,p,t)
    Rv, Ci = p
    Text = ext_disturbance(t)
    phih= ext_disturbance2(t)
    
    du[1] = 1/(Rv*Ci) .* (Text .- u[1]) .+ phih/Ci

end 

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

Text_1minuto= data[1:60]
 
Testerna = Text_1minuto;
phi1=0*ones(10)
phi2=0*ones(10)
phi3=0.5*ones(10)
phi4=1*ones(10)
phi5=0.5*ones(10)
phi6=0*ones(10)

phih=vcat(phi1,phi2,phi3,phi4,phi5,phi6);

tspan= (0.0f0, 60.0f00)
t= range(tspan[1], tspan[2], length=60)

disturbance = LinearInterpolation(Testerna,t);
disturbance2 = LinearInterpolation(phih,t);
u0= [4.0]

p= [0.05, 10.0]

function ext_disturbance(t)
return disturbance(t)
end

function ext_disturbance2(t)
return disturbance2(t)
end 
    

prob= ODEProblem(RC!, u0, tspan, p)
sol= Array(solve(prob, Tsit5(), saveat=t))
plot(t,sol[1,:], label = "reale",xlabel = "tempo [s]", ylabel = "Temperatura [°C]")

p = [0.05, 10]
rng= Random.default_rng()

dudt = Lux.Chain( x-> x,
Lux.Dense(1,50, tanh),
Lux.Dense(50,1))

p, st = Lux.setup(rng,dudt)
n_ode= NeuralODE(dudt, tspan, Tsit5(), saveat = t)

function predict_neuralode(p)
    Array(n_ode(u0,p,st)[1])
end 

function loss_neuralode(p)
    pred= predict_neuralode(p)
    loss = sum(abs2,  sol .- pred)
    return loss, pred
end 

callback = function (p, l, pred; doplot = false)
  println(l)
  # plot current prediction against data
  if doplot
    plt = plot(t, sol[1,:], label = "data")
    plot!(plt, t, pred[1,:], label = "prediction")
    display(plot(plt))
  end
  return false
end

pinit = ComponentArray(p)
callback(pinit, loss_neuralode(pinit)...)


adtype = Optimization.AutoZygote()

optf = Optimization.OptimizationFunction((x, p) -> loss_neuralode(x), adtype)
optprob = Optimization.OptimizationProblem(optf, pinit)


result_neuralode = Optimization.solve(optprob,
                                       Optimisers.Adam(0.05),
                                       callback = callback,
                                       maxiters = 400)
callback(result_neuralode.u, loss_neuralode(result_neuralode.u)...; doplot=true)


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


result_neuralode2 = Optimization.solve(optprob2,
                                     Optim.BFGS(initial_stepnorm=0.01),
                                        callback = callback,
                                        allow_f_increases = false)

#allow_f_increases impone all'ottimizzazione di stopparsi in automatico nel momento in cui è vicina al minimo

callback(result_neuralode2.u, loss_neuralode(result_neuralode2.u)...; doplot=true)
5

Did you solve the problem? using neural ode I always have poor learning results

6

Normally if you’re just using a pure neural ODE then it’s a bad idea. You want to build a model that has more predictive information in it unless you have tons of data. This of course is a completely different problem so if you have questions please ask in a new thread.