I have written a battery model in Julia. This model already runs without any issues.
I want to now do parameter estimation w.r.t to testing data using DiffEqFlux.jl . However, whenever I run the diffeqflux.scimltrain command, julia crashes with this message
The terminal process “julia.exe ‘-i’, ‘–banner=no’, ‘–project=C:\Users\madmin1.julia\environments\v1.6’, ‘c:\Users\madmin1.vscode\extensions\julialang.language-julia-1.3.34\scripts\terminalserver\terminalserver.jl’, ‘\.\pipe\vsc-jl-repl-76ed54af-dfa3-4f44-b40b-ceee81f8fccd’, ‘\.\pipe\vsc-jl-cr-2d720729-fad4-47fe-9903-02a83bb0f643’, ‘USE_REVISE=true’, ‘USE_PLOTPANE=true’, ‘USE_PROGRESS=true’, ‘DEBUG_MODE=false’” terminated with exit code: 3221226356.
I have already run everything else. Its the sciml line that’s causing the Julia to crash. I have also tried restarting the Julia language server and vscode too.
I am running Windows Server 2019 and have the latest Visual studio and Julia versions respectively.
‘’’
'''
using DifferentialEquations
using Plots
using DataFrames
using DiffEqFlux, Optim
using CSV
using Flux
plotly()
global const F = 96485.0 # Faraday's constant
global const R = 8.314462 # J.K⁻¹.mol⁻¹
T = 298 # K
function Δ!(A,N,D,R,Δr)
"""
Creates the 2nd order FDM matrix for solid phase concentration equation with Neumann boundary conditions
Comes from the following paper https://drive.google.com/file/d/1czf_P0tt8ViULMVcxTla_u-zaiGCsRE-/view?usp=sharing Page 41
"""
#Neumann BCs
A[1,1] = -D*(R[2] + Δr)/(R[2] *Δr^2)
A[1,2] = D*(R[2] + Δr)/(R[2] *Δr^2)
A[N,N-1] = D*(R[N-1] - Δr)/(R[N-1] *Δr^2)
A[N,N] = -D*(R[N-1] - Δr)/(R[N-1] *Δr^2)
for i in 2:N-1
A[i,i] = -2 * D/ (Δr^2)
A[i,i-1] = D*(R[i] - Δr)/(R[i] *Δr^2)
A[i,i+1] = D*(R[i] + Δr)/(R[i] *Δr^2)
end
end
function U₊_vs_c(θ)
"""
Returns Cathode equilibrium potential. Source: https://drive.google.com/file/d/16WE-LDKcFBVcLXrl-y-Y5UyfWY1u4GvE/view?usp=sharing
"""
return 4.4875 - 0.8090*θ - 0.0428*tanh(18.5138*(θ-0.5542)) - 17.732*tanh(15.789*(θ-0.3117)) + 17.5842*tanh(15.9308*(θ-0.3120))
end
function U₋_vs_c(θ)
"""
Returns Anode equilibrium potential. Source: https://drive.google.com/file/d/16WE-LDKcFBVcLXrl-y-Y5UyfWY1u4GvE/view?usp=sharing
"""
return 0.2482 - 0.0909*tanh(29.8538*(θ-0.1234)) - 0.04478*tanh(14.9159*(θ-0.2769)) - 0.0205*tanh(30.4444*(θ-0.6103)) + 1.9793*exp(-39.3631*θ)
end
function Butler_Volmer(c⁺ₛ , c⁻ₛ , p,T=298,I=6.5)
"""
Returns the cell voltage using Butler_Volmer equation. From paper: https://drive.google.com/file/d/1KwhI43EkDnFEZ2QGybs4o6bq7SOnPugx/view?usp=sharing
"""
R⁺ₛ ,R⁻ₛ , ε₊ , ε₋ , A , L₊ , L₋ , D₊ , D₋ , k⁺ , k⁻ , C⁰Lᵢ , cₘₐₓ⁺ , cₘₐₓ⁻ ,Rᶠ = p
α⁺ = α⁻ = 0.5
try
i₀⁺ = k⁺ * sqrt(C⁰Lᵢ * c⁺ₛ * (cₘₐₓ⁺ - c⁺ₛ ))
i₀⁻ = k⁻ * sqrt(C⁰Lᵢ * c⁻ₛ * (cₘₐₓ⁻ - c⁻ₛ ))
return (R*T/(α⁺*F)) * asinh(I*R⁺ₛ/(2*ε₊*A*i₀⁺)) - (R*T/(α⁻*F)) * asinh(I*R⁻ₛ/(2*ε₋*A*i₀⁻)) + U₊_vs_c(c⁺ₛ/cₘₐₓ⁺) - U₋_vs_c(c⁻ₛ/cₘₐₓ⁻) - (Rᶠ * I)
catch
return 0.0# Cell_parameters["Lower cutoff voltage"]
end
end
#Initial Conditions
Cₛ⁺ = ones(Float64,50) .* 17038.0 # calculated concentrations for Cathode
Cₛ⁻ = ones(Float64,50) .* 29866.0 # calculated concentrations for Anode
dCₜ = 0.0 .* ArrayPartition(Cₛ⁺,Cₛ⁻) #Intial fluxes are zero
Cₛ = ArrayPartition(Cₛ⁺,Cₛ⁻)
function SPM!(dCₜ,Cₛ,p,t)
Cell_R⁺ₛ ,Cell_R⁻ₛ , ε₊ , ε₋ , A , L₊ , L₋ , D₊ , D₋ , k⁺ , k⁻ , C⁰Lᵢ , cₘₐₓ⁺ , cₘₐₓ⁻ ,Rᶠ = p
I = 6.5
N₊ = 50
N₋ = 50
Δr₊ = Cell_R⁺ₛ/(N₊ - 1)
Δr₋ = Cell_R⁻ₛ/(N₋ -1)
Rs⁺ₛ = collect(0.0:Δr₊:Cell_R⁺ₛ)
Rs⁻ₛ = collect(0.0:Δr₋:Cell_R⁻ₛ)
Acs⁺ = zeros(Float64,N₊,N₊) #FDM matrix for positive solid phase concentration
Acs⁻ = zeros(Float64,N₋,N₋) #FDM matrix for negative solid phase concentration
Bcs⁺ = zeros(Float64,N₊,1)
Bcs⁺[N₊] = (Rs⁺ₛ[N₊ - 1] + Δr₊ )*Cell_R⁺ₛ/(Rs⁺ₛ[N₊-1] * Δr₊ * 3 * ε₊ * F * A * L₊) #Neumann BC source term for positive
Bcs⁻ = zeros(Float64,N₋,1)
Bcs⁻[N₋] = -(Rs⁻ₛ[N₋ - 1] + Δr₋ )*Cell_R⁻ₛ/(Rs⁻ₛ[N₋ - 1] * Δr₋ * 3 * ε₋ * F * A * L₋) #Neumann BC source term for negative
Δ!(Acs⁺,N₊,D₊,Rs⁺ₛ,Δr₊)
Δ!(Acs⁻,N₋,D₋,Rs⁻ₛ,Δr₋)
dCₜ[1:N₊] = Acs⁺ * Cₛ[1:N₊] + Bcs⁺ .* I
dCₜ[N₊+1:N₋+N₊] = Acs⁻ * Cₛ[N₊+1:N₋+N₊] + Bcs⁻ .* I
end
tspan = (0.0,2300)
p = [5.22e-6,5.86e-6,0.75,0.665,0.1027,7.56e-5,8.52e-5,4.0e-15,3.3e-14,3.42e-6,6.48e-7,1000,63104.0,33133.0,0.1]
SPM_model = ODEProblem(SPM!,Cₛ,tspan,p)
sol = solve(SPM_model,alg_hints = [:stiff])
ts_BV = collect(sol.t[1]:1:sol.t[end])
V1=[]
for t in ts_BV
push!(V1,Butler_Volmer(sol(t)[50],sol(t)[100],p))
end
plot(ts_BV,V1,label="Initial SPM Prediction")
df = DataFrame(CSV.File("D:\\Shwetank files\\Google Drive\\Life_prediction_data_sets\\CSV_data\\SNL_18650_NCA_25C_0-100_0.5-2C_b_timeseries.csv"))
df = df[!,["Test_Time (s)","Cycle_Index","Current (A)" , "Voltage (V)"]]
dropmissing(df)
filter!(row -> row."Current (A)" < 0.0,df)
Cycle1 = filter(row -> row."Cycle_Index"==6.0,df)
Cycle1[!,"Current (A)"] = Cycle1[!,"Current (A)"].*-1
Cycle1[!,"Test_Time (s)"] = Cycle1[!,"Test_Time (s)"] .- Cycle1[!,"Test_Time (s)"][1]
plot!(Cycle1[!,"Test_Time (s)"],Cycle1[!,"Voltage (V)"],label = "SNL_1C")
function loss(p)
tmp_prob = remake(SPM_model,p=p)
tmp_sol = solve(tmp_prob,alg_hints = [:stiff])
V1 = []
for t in Cycle1[!,"Test_Time (s)"]
push!(V1,Butler_Volmer(tmp_sol(t)[50],tmp_sol(t)[100],p))
end
l = sum(abs2,V1 - Cycle1[!,"Voltage (V)"]) |> sqrt
return l, tmp_sol
end
pinit = p
cb= function (p,l,pred)
println("Loss: $l")
return false
end
res = DiffEqFlux.sciml_train(loss,pinit,ADAM(0.01),cb=cb,maxiters=10)
julia> versioninfo()
Julia Version 1.6.2
Commit 1b93d53fc4 (2021-07-14 15:36 UTC)
Platform Info:
OS: Windows (x86_64-w64-mingw32)
CPU: Intel(R) Xeon(R) Gold 5220 CPU @ 2.20GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-11.0.1 (ORCJIT, cascadelake)
Environment:
JULIA_EDITOR = code
JULIA_NUM_THREADS =