Hi all,
I’m following these Universal ODE examples, trying to get one of my own systems to work. Having a good deal of issues getting things to solve & I hope someone can see the error I’m making.
Problem is that I see no errors turning up, some of the function calls are just spooling forever. I’ve turned on status bars on the ODE solver calls and I either get
ODE 0%| | ETA: N/A
or stuck somewhere around 25% with the ETA just continually climbing up, never going down.
My system:
Where U_{1,2} are the (arbitrary) portions of my system I’d like to try and recover. They are m_p P and c_V V^2 respectively in truth.
The tspan for my problem is large compared to the L-V example. The system runs on a time unit of days and the dynamics I’m interested in cover 30 years (10,950 days). A normal solve with Tsit5() is practically instantaneous over that range.
I’ve set up the training set to cover 3 years, with 30 data points over that span.
Current code:
using OrdinaryDiffEq
using ModelingToolkit
using DataDrivenDiffEq
using LinearAlgebra, DiffEqSensitivity, Optim
using DiffEqFlux, Flux
using JLD2
println("Generate data")
function lake!(du, u, p, t)
B, P, V = u
N, ib, ip, r, H₁, H₂, H₃, H₄, cb, cp, prmax, ce, mp, rv, cv, mv = p
FR = B^2 / (B^2 + H₄^2)
du[1] =
dB = ib + r * (N / (N + H₁)) * B - cb * B^2 - prmax * FR * P
du[2] =
dP = ip + ce * prmax * FR * P * (V / (V + H₂)) - mp * P - cp * P^2
du[3] =
dV = rv * V - cv * V^2 - mv * (V * B^2 / (H₃^2 + B^2))
end
# Define the experiment
datasize = 30
tspan = (0.0f0, 1095.0f0) # Full data would be 10950: 30 years in days. Currently 3 years with 30 training points
tsteps = range(tspan[1], tspan[2], length = datasize)
u0 = Float32[60.036, 0.738, 11.654] # X1
p_ = Float32[
0.7,
2e-5,
2e-5,
7.5e-3,
0.5,
20,
20,
15,
7.5e-5,
2.75e-4,
5e-2,
0.14,
2.25e-3,
7e-3,
6e-5,
7e-3,
]
prob = ODEProblem(lake!, u0, tspan, p_)
solution = solve(prob, Vern7(), abstol = 1e-12, reltol = 1e-12, maxiters = 1e7, saveat = tsteps)
# Ideal data
X = Array(solution)
# Add noise to the ideal data
Xₙ = X + Float32(1e-3) * randn(eltype(X), size(X))
# Training implementation
# Input layer (BPV), hidden layer, ouptput layer (u1, u2)
ann_L = FastChain(FastDense(3, 32, tanh), FastDense(32, 32, tanh), FastDense(32, 2))
p = initial_params(ann_L)
function dudt_(u, p, t)
B, P, V = u
N, ib, ip, r, H₁, H₂, H₃, H₄, cb, cp, prmax, ce, mp, rv, cv, mv = p_
# z should be u1, u2 correct?
z = ann_L(u, p)
FR = B^2 / (B^2 + H₄^2)
[ib + r * (N / (N + H₁)) * B - cb * B^2 - prmax * FR * P,
ip + ce * prmax * FR * P * (V / (V + H₂)) - z[1] - cp * P^2,
rv * V - z[2] - mv * (V * B^2 / (H₃^2 + B^2))]
end
prob_nn = ODEProblem(dudt_, u0, tspan, p)
# contrete_solve has been depreciated apparently, although line below fails if I just use solve at present
# Depending on initialization, this line is the one that gets a portion of the way then just stops.
# Other times it solves almost instantly.
sol_nn = concrete_solve(prob_nn, Tsit5(), u0, p, saveat = solution.t)
function predict(θ)
Array(concrete_solve(
prob_nn,
Vern7(),
u0,
θ,
saveat = solution.t,
abstol = 1e-6,
reltol = 1e-6,
maxiters = 1e7,
sensealg = InterpolatingAdjoint(autojacvec = ReverseDiffVJP()),
))
end
function loss(θ)
pred = predict(θ)
sum(abs2, Xₙ .- pred), pred
end
# Test - this works without issue so long as sol_nn does.
loss(p)
losses = []
callback(θ, l, pred) = begin
push!(losses, l)
if !isempty(losses)
println("Current loss after $(length(losses)) iterations: $(losses[end])")
end
false
end
# First train with ADAM for better convergence
res1 = DiffEqFlux.sciml_train(loss, p, ADAM(0.01), cb=callback, maxiters=200)
# Rest of the script is irrelevant for us.
And here is where I run into issues. res1 never starts: the status bar just sits at Training 0%| | ETA: N/A which gives me little to go on in terms of troubleshooting - there’s a lot of packages I’m unfamiliar with and the problem could be in any of them, or more likely, my set-up.
The system doesn’t seem to be stiff. I’ve tried with AutoTsit5(Rosenbrock23()) and AutoVern7(Rodas5()) just in case - no difference.
I’ve attempted to solve this on the gpu as well, by
u0 |> gpu
p |> gpu
p_ |> gpu
res1 = DiffEqFlux.sciml_train(loss, p, ADAM(0.01), cb=callback, maxiters=200)
but that also doesn’t do much of anything.
Can you see the issue at all?
