Error during Optimization using Neural ODEs

Marco_Nesta · October 17, 2023, 4:20pm

using DifferentialEquations, Plots, Flux,Optim, DiffEqFlux, DataInterpolations,Random, ComponentArrays, Lux
using Optimization, OptimizationOptimisers, OptimizationOptimJL,OptimizationNLopt
rng = Random.default_rng()
using CSV
using DataFrames
using Plots
using Flux
using Optim: LineSearches, BFGS
using Statistics: mean
using DiffEqFlux, Optimization, OptimizationOptimJL,Plots
using ComponentArrays, Lux, DiffEqFlux, Optimization, OptimizationPolyalgorithms, DifferentialEquations, Plots
using DiffEqFlux: group_ranges
using WeightInitializers: truncated_normal
# Load the data
df = CSV.read("es_01b_real_data_.csv", DataFrame)
df=repeat(df, outer=1)

registered_gas_flow = df[1:10:4001, :2]
Registered_Temperature = df[1:10:4001, :3]
tsteps= df[1:10:4001, :1]

gas_flow = LinearInterpolation(registered_gas_flow,tsteps);
function ext_flow(tsteps)
    return gas_flow(tsteps)
end
function water_temp(tsteps)
    return Temperature_h20(tsteps)
end

gas_flow = LinearInterpolation(registered_gas_flow,tsteps);
Temperature_h20=LinearInterpolation(Registered_Temperature,tsteps);


function RC!(du,u,p,t)
    x = u[1]
    A,B,C,D = p
    P= ext_flow(t)
    du[1]=dx=(B/A)* P - (C/A) + (D * (20 -x))/A
    end

u0= [20.0]

tspan= (0.0f0,400.0f0)
global p = [66.896,50e6,100.0,0.2]; # p, M, e

prob= ODEProblem(RC!, u0, tspan, p)
ode_data=Array(solve(prob,saveat=tsteps))
#normalize the data between 0 and 1
ode_data=ode_data./maximum(ode_data)
ode_data=ode_data.-minimum(ode_data)

plot(ode_data[1,:], label="ODE")

const nn = Lux.Chain(Base.Fix1(broadcast, cos),
    Lux.Dense(1 => 32, cos; init_weight=truncated_normal(; std=1e-4)),
    Lux.Dense(32 => 32, cos; init_weight=truncated_normal(; std=1e-4)),
    Lux.Dense(32 => 1; init_weight=truncated_normal(; std=1e-4)))
ps, st = Lux.setup(MersenneTwister(), nn)

params = ComponentArray{Float64}(ps)

function ODE_model(u, nn_params, t)
    P= ext_flow(t)
    A,B,C = p
    y =(B/A).* P .- (100/66.896) .+ first(nn([first(u)], nn_params, st))
end

prob_nn = ODEProblem(ODE_model, u0, tspan, params)
soln_nn = Array(solve(prob_nn,saveat=1.0))
soln_nn=soln_nn./maximum(soln_nn)
soln_nn=soln_nn.-minimum(soln_nn)

plot(soln_nn[1,:], label="NN")
plot!(ode_data[1,:], label="ODE")

function loss(θ)
    pred = Array(solve(prob_nn,saveat=tsteps))
    loss = sum(abs2, ode_data[1,:] .- soln_nn[1,:])
    return loss, pred
end

loss(params)

losses = Float64[]

function callback(θ, l, pred)
    push!(losses, l)
    println("Training || Iteration: $(length(losses)) || Loss: $(l)")
    return false
end

adtype = Optimization.AutoZygote()
optf = Optimization.OptimizationFunction((x, p) -> loss(x), adtype)
optprob = Optimization.OptimizationProblem(optf, params)
res = Optimization.solve(optprob,
    BFGS(; initial_stepnorm=0.01, linesearch=LineSearches.BackTracking());
    callback, maxiters=1000)

ERROR: MethodError: Cannot `convert` an object of type Nothing to an object of type Float64

Closest candidates are:
  convert(::Type{T}, ::Gray24) where T<:Real
   @ ColorTypes C:\Users\marco\.julia\packages\ColorTypes\1dGw6\src\conversions.jl:114
  convert(::Type{T}, ::Gray) where T<:Real
   @ ColorTypes C:\Users\marco\.julia\packages\ColorTypes\1dGw6\src\conversions.jl:113
  convert(::Type{T}, ::LLVM.GenericValue, ::LLVM.LLVMType) where T<:AbstractFloat
   @ LLVM C:\Users\marco\.julia\packages\LLVM\HykgZ\src\execution.jl:39
  ...

Stacktrace:
  [1] fill!(A::ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}}}, x::Nothing)
    @ Base .\multidimensional.jl:1084
  [2] copyto!
    @ .\broadcast.jl:934 [inlined]
  [3] materialize!
    @ .\broadcast.jl:884 [inlined]
  [4] materialize!(dest::ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}}}, bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0}, Nothing, typeof(identity), Tuple{Base.RefValue{Nothing}}})
    @ Base.Broadcast .\broadcast.jl:881
  [5] (::Optimization.var"#293#302"{Optimization.var"#292#301"{OptimizationFunction{true, Optimization.AutoZygote, var"#11#12", Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Optimization.ReInitCache{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}}}, SciMLBase.NullParameters}}})(::ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}}}, ::ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}}})
    @ Optimization C:\Users\marco\.julia\packages\Optimization\vFala\src\function\zygote.jl:74
  [6] (::OptimizationOptimJL.var"#10#16"{OptimizationOptimJL.OptimJLOptimizationCache{OptimizationFunction{false, Optimization.AutoZygote, OptimizationFunction{true, Optimization.AutoZygote, var"#11#12", Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Optimization.var"#293#302"{Optimization.var"#292#301"{OptimizationFunction{true, Optimization.AutoZygote, var"#11#12", Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Optimization.ReInitCache{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}}}, SciMLBase.NullParameters}}}, Optimization.var"#296#305"{Optimization.var"#292#301"{OptimizationFunction{true, Optimization.AutoZygote, var"#11#12", Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Optimization.ReInitCache{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}}}, SciMLBase.NullParameters}}}, Optimization.var"#300#309", Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Optimization.ReInitCache{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}}}, SciMLBase.NullParameters}, Nothing, Nothing, Nothing, Nothing, Nothing, BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Flat}, Base.Iterators.Cycle{Tuple{Optimization.NullData}}, Bool, typeof(callback)}, OptimizationOptimJL.var"#9#15"{OptimizationOptimJL.OptimJLOptimizationCache{OptimizationFunction{false, Optimization.AutoZygote, OptimizationFunction{true, Optimization.AutoZygote, var"#11#12", Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Optimization.var"#293#302"{Optimization.var"#292#301"{OptimizationFunction{true, Optimization.AutoZygote, var"#11#12", Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Optimization.ReInitCache{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}}}, SciMLBase.NullParameters}}}, Optimization.var"#296#305"{Optimization.var"#292#301"{OptimizationFunction{true, Optimization.AutoZygote, var"#11#12", Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Optimization.ReInitCache{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}}}, SciMLBase.NullParameters}}}, Optimization.var"#300#309", Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Optimization.ReInitCache{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}}}, SciMLBase.NullParameters}, Nothing, Nothing, Nothing, Nothing, Nothing, BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Flat}, Base.Iterators.Cycle{Tuple{Optimization.NullData}}, Bool, typeof(callback)}}})(G::ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}}}, θ::ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}}}) 
    @ OptimizationOptimJL C:\Users\marco\.julia\packages\OptimizationOptimJL\uRfW9\src\OptimizationOptimJL.jl:200
  [7] value_gradient!!(obj::TwiceDifferentiable{Float64, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}}}, ComponentMatrix{Float64, Matrix{Float64}, Tuple{Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}, Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}}}, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}}}}, x::ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}}})
    @ NLSolversBase C:\Users\marco\.julia\packages\NLSolversBase\kavn7\src\interface.jl:82
  [8] initial_state(method::BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Flat}, options::Optim.Options{Float64, OptimizationOptimJL.var"#_cb#14"{OptimizationOptimJL.OptimJLOptimizationCache{OptimizationFunction{false, Optimization.AutoZygote, OptimizationFunction{true, Optimization.AutoZygote, var"#11#12", Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Optimization.var"#293#302"{Optimization.var"#292#301"{OptimizationFunction{true, Optimization.AutoZygote, var"#11#12", Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Optimization.ReInitCache{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}}}, SciMLBase.NullParameters}}}, Optimization.var"#296#305"{Optimization.var"#292#301"{OptimizationFunction{true, Optimization.AutoZygote, var"#11#12", Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Optimization.ReInitCache{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}}}, SciMLBase.NullParameters}}}, Optimization.var"#300#309", Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Optimization.ReInitCache{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}}}, SciMLBase.NullParameters}, Nothing, Nothing, Nothing, Nothing, Nothing, BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Flat}, Base.Iterators.Cycle{Tuple{Optimization.NullData}}, Bool, typeof(callback)}}}, d::TwiceDifferentiable{Float64, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}}}, ComponentMatrix{Float64, Matrix{Float64}, Tuple{Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}, Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}}}, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}}}}, initial_x::ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}}})
    @ Optim C:\Users\marco\.julia\packages\Optim\V8ZEC\src\multivariate\solvers\first_order\bfgs.jl:94
  [9] optimize(d::TwiceDifferentiable{Float64, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}}}, ComponentMatrix{Float64, Matrix{Float64}, Tuple{Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}, Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}}}, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}}}}, initial_x::ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}}}, method::BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Flat}, options::Optim.Options{Float64, OptimizationOptimJL.var"#_cb#14"{OptimizationOptimJL.OptimJLOptimizationCache{OptimizationFunction{false, Optimization.AutoZygote, OptimizationFunction{true, Optimization.AutoZygote, var"#11#12", Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Optimization.var"#293#302"{Optimization.var"#292#301"{OptimizationFunction{true, Optimization.AutoZygote, var"#11#12", Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Optimization.ReInitCache{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}}}, SciMLBase.NullParameters}}}, Optimization.var"#296#305"{Optimization.var"#292#301"{OptimizationFunction{true, Optimization.AutoZygote, var"#11#12", Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Optimization.ReInitCache{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}}}, SciMLBase.NullParameters}}}, Optimization.var"#300#309", Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Optimization.ReInitCache{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}}}, SciMLBase.NullParameters}, Nothing, Nothing, Nothing, Nothing, Nothing, BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Flat}, Base.Iterators.Cycle{Tuple{Optimization.NullData}}, Bool, typeof(callback)}}})
    @ Optim C:\Users\marco\.julia\packages\Optim\V8ZEC\src\multivariate\optimize\optimize.jl:36
 [10] __solve(cache::OptimizationOptimJL.OptimJLOptimizationCache{OptimizationFunction{false, Optimization.AutoZygote, OptimizationFunction{true, Optimization.AutoZygote, var"#11#12", Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Optimization.var"#293#302"{Optimization.var"#292#301"{OptimizationFunction{true, Optimization.AutoZygote, var"#11#12", Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Optimization.ReInitCache{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(layer_1 = 1:0, layer_2 = ViewAxis(1:64, Axis(weight = ViewAxis(1:32, ShapedAxis((32, 1), NamedTuple())), bias = ViewAxis(33:64, ShapedAxis((32, 1), NamedTuple())))), layer_3 = ViewAxis(65:1120, Axis(weight = ViewAxis(1:1024, ShapedAxis((32, 32), NamedTuple())), bias = ViewAxis(1025:1056, ShapedAxis((32, 1), NamedTuple())))), layer_4 = ViewAxis(1121:1153, Axis(weight = ViewAxis(1:32, ShapedAxis((1, 32), NamedTuple())), bias = ViewAxis(33:33, ShapedAxis((1, 1), NamedTuple())))))}}}, SciMLBase.NullParameters}}}, Optimization.var"#296#305"{Optimization.var"#292#301"{OptimizationFunction{true, Optimization.AutoZygote, var"#11#12", Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing,

contradict · October 17, 2023, 4:43pm

I’m not sure this is your problem, but it appears that loss is independent of θ. You probably want something like

function loss(θ)
    prob = remake(prob_nn; p=θ)
    pred = Array(solve(prob,saveat=tsteps))
    loss = sum(abs2, ode_data[1,:] .- pred[1,:])
    return loss, pred
end

Marco_Nesta · October 18, 2023, 2:39pm

Thank you. the code seems to work but I continue to have low prediction accuracy. Is there something I am doing wrong? I want to use neural odes to learn the dynamics of the model I define. The model has 4 parameters and two exogenous variables

using DifferentialEquations, Plots, Flux,Optim, DiffEqFlux, DataInterpolations,Random, ComponentArrays, Lux
using Optimization, OptimizationOptimisers, OptimizationOptimJL,OptimizationNLopt
rng = Random.default_rng()
using CSV
using DataFrames
using Plots
using Flux
using Optim: LineSearches, BFGS
using Statistics: mean,sqrt!,abs2
using DiffEqFlux, Optimization, OptimizationOptimJL,Plots
using ComponentArrays, Lux, DiffEqFlux, Optimization, OptimizationPolyalgorithms, DifferentialEquations, Plots
using DiffEqFlux: group_ranges
using WeightInitializers: truncated_normal
# Load the data
df = CSV.read("es_02_real_data_.csv", DataFrame)
df=repeat(df, outer=1)
end_param= 18001
tsteps= df[1:10:end_param, :1]
registered_gas_flow = df[1:10:end_param, :2]

T_ext= df[1:10:end_param, :3]
T_h20= df[1:10:end_param, :4]   

function ext_Temp(tsteps)
    return Temperature_ext(tsteps)
end
    

function ext_flow(tsteps)
    return gas_flow(tsteps)
end
function water_temp(tsteps)
    return Temperature_h20(tsteps)
end

gas_flow = LinearInterpolation(registered_gas_flow,tsteps);
Temperature_h20=LinearInterpolation(T_h20,tsteps);
Temperature_ext=LinearInterpolation(T_ext,tsteps);


#create a 1800 time vector
function RC!(du,u,p,t)
    x = u[1]
    A,B,C,D = p
    P= ext_flow(t)
    Testerna=ext_Temp(t)
    du[1]=dx=(B/A)* P - (C/A) + (D * (Testerna-x))/A
end

u0= [20.0]

tspan= (0.0f0,1800.0f0)
global p = [66.896,50e6,100.0,3.3]; # p, M, e

gas_flow = LinearInterpolation(registered_gas_flow,tsteps);
Temperature_h20=LinearInterpolation(T_h20,tsteps);
Temperature_ext=LinearInterpolation(T_ext,tsteps);

prob= ODEProblem(RC!, u0, tspan, p)

u0= [20.0]

ode_data=Array(solve(prob,saveat=tsteps))


test_data_gas = df[18001:10:36001, :2]
tsteps_test= df[18001:10:36001, :1]
test_data_text=df[18001:10:36001, :3]

gas_flow_test = LinearInterpolation(test_data_gas,tsteps_test);
Temperature_ext_test=LinearInterpolation(test_data_text,tsteps_test);
function ext_flow_test(tsteps_test)
    return gas_flow_test(tsteps_test)
end

function ext_Temp_test(tsteps_test)
    return Temperature_ext_test(tsteps_test)
end

#create a 3600 time vector
function RC_test(du,u,p,t)
    x = u[1]
    A,B,C,D = p
    P= ext_flow_test(t)
    Text= ext_Temp_test(t)
    du[1]=dx=1/A*(B* P - C + D * (Text -x))
    end


tspan_test= (tspan[2],3600.0f0)
u1=[ode_data[1,end]]


prob_test= ODEProblem(RC_test, u1, tspan_test, p)
ode_data_test=Array(solve(prob_test,saveat=tsteps_test))

plot(tsteps,ode_data[1,:], label="ode_data")
plot!(tsteps_test,ode_data_test[1,:], label="ode_data_test")
const nn = Lux.Chain(Base.Fix1(broadcast, sin),
    Lux.Dense(1 => 32, sin; init_weight=truncated_normal(; std=1e-4)),
    Lux.Dense(32 => 32, sin; init_weight=truncated_normal(; std=1e-4)),
    Lux.Dense(32 => 1; init_weight=truncated_normal(; std=1e-4)))
ps, st = Lux.setup(MersenneTwister(), nn)

const params = ComponentArray{Float64}(ps)

function ODE_model(u, nn_params, t)
    P= ext_flow(t)
    Testerna=ext_Temp(t)
    A,B,C,D = p
    y =1/A.*(B.* P .- C .+ D.*(first(nn([first(u)], nn_params, st))))
end

prob_nn = ODEProblem(ODE_model, u0, tspan, params)
soln_nn = Array(solve(prob_nn,saveat=tsteps))

plot(soln_nn[1,:], label="NN")
plot!(ode_data[1,:], label="ODE")

function loss(θ)
    pred = Array(solve(prob_nn, Tsit5(); u0, p=θ, saveat=tsteps))
    loss = sqrt(mean(abs2.(ode_data[1,:] .- pred[1,:])))
    return loss, pred
end

loss(params)

const losses = Float64[]

function callback(θ, l, pred)
    push!(losses, l)
    println("Training || Iteration: $(length(losses)) || Loss: $(l)")
    return false
end

adtype = Optimization.AutoZygote()
optf = Optimization.OptimizationFunction((x, p) -> loss(x), adtype)
optprob = Optimization.OptimizationProblem(optf, params)
res = Optimization.solve(optprob,Optimisers.Adam(); callback, maxiters=100)
optprob2 = remake(optprob, u0 = res.u)
res2 = Optimization.solve(optprob2,
    BFGS(; initial_stepnorm=0.01, linesearch=LineSearches.BackTracking());
    callback, maxiters=1000)

plot(losses)
display(plot!(xlabel="Iteration", ylabel="Loss"))

prob_nn = ODEProblem(ODE_model, u0, tspan, res2.u)
soln_nn =Array(solve(prob_nn, Tsit5(); u0, p=res2.u, saveat=tsteps))

plot(soln_nn[1,:], label="NN trained")
display(plot!(ode_data[1,:], label="ODE"))

ChrisRackauckas · October 18, 2023, 10:12pm

Reduce the tolerances on the ODE solver and try things like multiple shooting?

Marco_Nesta · October 19, 2023, 7:57am

Thanks Chris, I’ll try. what do you think of the code? does it seem like a right approach? sorry but I’m new to this approach

ChrisRackauckas · October 19, 2023, 9:16am

It’s the right approach, but tolerances matter.

Marco_Nesta · October 19, 2023, 10:34am

Can normalizing the data help?

ChrisRackauckas · October 19, 2023, 10:37am

yes it definitely can

Marco_Nesta · October 19, 2023, 10:46am

my doubt was that by normalizing the parameters, the dynamics were not learned in the right way

Marco_Nesta · October 19, 2023, 10:48am

nntrained
this is the solution adding the tollerances in ode solvers

Marco_Nesta · October 19, 2023, 10:49am

I think the error is due to the fact that most of the variation is given by Text-Th20 which is what my network has to learn, and which it fails to do

ChrisRackauckas · October 19, 2023, 10:50am

And then multiple shooting? Show the data points on the graph?

Marco_Nesta · October 19, 2023, 12:53pm

poits

ChrisRackauckas · October 19, 2023, 1:29pm

and what happens with multiple shooting?

Marco_Nesta · October 19, 2023, 1:42pm

I am now trying to decrease the density of the data so as to simplify the problem

Marco_Nesta · October 19, 2023, 1:50pm

initial nn
this is the first initialization of the network before optimization

ChrisRackauckas · October 19, 2023, 2:17pm

What are the extra start points you used for the multiple shooting?

Marco_Nesta · October 19, 2023, 2:22pm

now I am not using multiple shooting yet, I am using this link as a reference

ChrisRackauckas · October 19, 2023, 4:29pm

For this type of problem, multiple shooting or PEM is pretty required to get a good fit.

baggepinnen · October 19, 2023, 5:29pm

Is this data simulated? A simple linear model of order 2 fits perfectly

using DelimitedFiles
using ControlSystemIdentification
using Plots

data = readdlm("/tmp/es_02_real_data_.csv", ',')

t = data[2:end, 1] .|> Float64
u = data[2:end, 2] .|> Float64
u2 = data[2:end, 3] .|> Float64
y = data[2:end, 4] .|> Float64

Ts = t[2] - t[1]

d = iddata(y', [u u2]', Ts)
model = subspaceid(d, 2, scaleU=true, verbose=true)
simplot(model, d)

Or, if you’d like, train on first half and predict on second half of data

model = subspaceid(d[1:end÷2], 2, scaleU=true, verbose=true)
simplot(model, d[end÷2:end])

Maybe you aren’t supposed to use both input signals in the model?