using OrdinaryDiffEq , SciMLSensitivity ,Optimization, OptimizationOptimisers,OptimizationOptimJL, LineSearches
using Statistics
using StableRNGs, Lux, Zygote , Plots , ComponentArrays
rng = StableRNG(11)
# Generating training data
function actualODE!(du,u,p,t,T∞,I)
Cbat = 5*3600
du[1] = -I/Cbat
C₁ = -0.00153 # Unit is s-1
C₂ = 0.020306 # Unit is K/J
R0 = 0.03 # Resistance set a 30mohm
Qgen =(I^2)*R0
du[2] = (C₁*(u[2]-T∞)) + (C₂*Qgen)
end
t1 = collect(0:1:3400)
T∞1,I1 = 298.15,5
actualODE1!(du,u,p,t) = actualODE!(du,u,p,t,T∞1,I1)
prob = ODEProblem(actualODE1!,[1.0,T∞1],(t1[1],t1[end]))
solution = solve(prob,Tsit5(),saveat = t1)
X = Array(solution)
T1 = X[2,:]
# Plotting the results
plot(solution[2,:],color = :red,label = ["True Data" nothing])
# Defining the neural network
const U = Lux.Chain(Lux.Dense(3,20,tanh),Lux.Dense(20,20,tanh),Lux.Dense(20,1))
_para,st = Lux.setup(rng,U)
const _st = st
function NODE_model!(du,u,p,t,T∞,I)
Cbat = 5*3600
du[1] = -I/Cbat
C₁ = -0.00153
C₂ = 0.020306
G = I*(U([u[1],u[2],I],p,_st)[1][1])
du[2] = (C₁*(u[2]-T∞)) + (C₂*G)
end
NODE_model1!(du,u,p,t) = NODE_model!(du,u,p,t,T∞1,I1)
prob1 = ODEProblem(NODE_model1!,[1.0,T∞1],(t1[1],t1[end]),_para)
function loss(θ)
_prob1 = remake(prob1,p=θ)
_sol = Array(solve(_prob1,Tsit5(),saveat = t1))
loss1 = mean(abs2,T1.-_sol[2,:])
return loss1
end
losses = Float64[]
callback = function(state,l)
push!(losses,l)
println("RMSE Loss at iteration $(length(losses)) is $sqrt(l)")
return false
end
adtype = Optimization.AutoZygote()
optf = Optimization.OptimizationFunction((x,p) -> loss(x),adtype)
optprob = Optimization.OptimizationProblem(optf,ComponentVector{Float64}(_para))
res1 = Optimization.solve(optprob, OptimizationOptimisers.Adam(),callback = callback,maxiters = 500)
Segfaults on v1.11: