Julia crashes during parameter estimation for a system of differential equations

[shwetankverma28]

I have asked this question in Julia discourse but didn't get any solution:
https://discourse.julialang.org/t/julia-crashes-when-doing-parameter-estimation-for-a-system-of-differential-equations-in-vscode/67859

Here's the same post:
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.

Here's the full code. I will also attach the CSV file with test data in case anyone wants to replicate it:

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
    
    """
    #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. 
    """
    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. 
    """
    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. 
    """
    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("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)

I have tried running this directly in Repl but it closes before I can read anything. I also ran this on a different PC and the same thing happened.
I ran the loss function separately and it works fine. I can't really figure what's going wrong.

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 = 60

SNL_18650_NCA_25C_0-100_0.5-2C_b_timeseries.csv

[shwetankverma28]

I ran this on Atom/Juno and this is what I got:

Please submit a bug report with steps to reproduce this fault, and any error messages that follow (in their entirety). Thanks.
Exception: EXCEPTION_ACCESS_VIOLATION at 0x7ff919a74408 -- RtlSetUserValueHeap at C:\Windows\SYSTEM32\ntdll.dll (unknown line)
in expression starting at C:\Users\verma\OneDrive\Desktop\Modelling\SPM.jl:140
RtlSetUserValueHeap at C:\Windows\SYSTEM32\ntdll.dll (unknown line)
RtlAllocateHeap at C:\Windows\SYSTEM32\ntdll.dll (unknown line)
RtlAllocateHeap at C:\Windows\SYSTEM32\ntdll.dll (unknown line)
malloc at C:\Windows\System32\msvcrt.dll (unknown line)

Please submit a bug report with steps to reproduce this fault, and any error messages that follow (in their entirety). Thanks.
Exception: EXCEPTION_ACCESS_VIOLATION at 0x7ff919a9e3cb --

julia> versioninfo()
Julia Version 1.6.2
Commit 1b93d53fc4 (2021-07-14 15:36 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)      
  CPU: AMD Ryzen 5 3500 6-Core Processor
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-11.0.1 (ORCJIT, znver2)
Environment:
  JULIA_EDITOR = "C:\Users\verma\AppData\Local\atom\app-1.58.0\atom.exe"  -a
  JULIA_NUM_THREADS = 6

I ran sfc /scannow and windows said that it fixed some corrupted files. I ran Atom/Juno again but this time it closed without giving message

[ChrisRackauckas]

Can you isolate it?

[shwetankverma28]

I am not sure how I would do that. I tried changing the loss function like this:

function loss(p)
    tmp_prob = remake(SPM_model,p=p)
    tmp_sol = solve(tmp_prob,alg_hints = [:stiff])
    test=ones(Float64,2000)
    V1 = []
    for t in 1:2000
        push!(V1,Butler_Volmer(tmp_sol(t)[50],tmp_sol(t)[100],p))
    end
    l = sum(abs2,V1 - test) |> sqrt
    return l, tmp_sol
end

And it still crashed. I thought it was reading from CSV that was causing the issue and but it still crashed.
I also tried just returning a constant loss value:

function loss(p)
    tmp_prob = remake(SPM_model,p=p)
    tmp_sol = solve(tmp_prob,alg_hints = [:stiff])
    test=ones(Float64,2000)
    V1 = []
    for t in 1:2000
        push!(V1,Butler_Volmer(tmp_sol(t)[50],tmp_sol(t)[100],p))
    end
    l = 1 #sum(abs2,V1 - test) |> sqrt
    return l, tmp_sol
end

I even tried this loss function:

function loss(p)
    tmp_prob = remake(SPM_model,p=p)
    tmp_sol = solve(tmp_prob,alg_hints = [:stiff])
 
    l = 1 #sum(abs2,V1 - test) |> sqrt
    return l, tmp_sol
end

And it still crashed. I am completely clueless here.

I should add that the original loss function works fine when I run it separately

[ChrisRackauckas]

Just keep making it simpler. I am sure not all of the details in the ODE are necessary? You can't delete a single line without it fixing the bug? Can't remove any constants? Etc. The way to figure out what's going on is to just iteratively make it simpler line by line.

[shwetankverma28]

I think I have isolated the issue to these four lines:

Δr₊ = Cell_R⁺ₛ/(N₊ - 1)
Δr₋ = Cell_R⁻ₛ/(N₋ -1)
Rs⁺ₛ = collect(0.0:Δr₊:Cell_R⁺ₛ)
Rs⁻ₛ = collect(0.0:Δr₋:Cell_R⁻ₛ)

Even if I comment out everything that uses these two arrays, it still crashes.

I am trying to estimate these parameters:
Cell_R⁺ₛ ,Cell_R⁻ₛ , ε₊ , ε₋ , A , L₊ , L₋ , D₊ , D₋ , k⁺ , k⁻ , C⁰Lᵢ , cₘₐₓ⁺ , cₘₐₓ⁻ ,Rᶠ = p
It's the first two parameters that are causing all the issues. Because if I don't estimate them and just hold them constant like this:

Cell_R⁺ₛ = 5.22e-6
Cell_R⁻ₛ = 5.86e-6
ε₊ , ε₋ , A , L₊ , L₋ , D₊ , D₋ , k⁺ , k⁻ , C⁰Lᵢ , cₘₐₓ⁺ , cₘₐₓ⁻ ,Rᶠ  = p

Then it doesn't crash. It does give me a warning though:
Warning: AD methods failed, using numerical differentiation. To debug, try ForwardDiff.gradient(loss, θ) or Zygote.gradient(loss, θ)
But it doesn't crash.

[shwetankverma28]

I have fixed the issue!
I replaced these two lines in ODE function:

Rs⁺ₛ = collect(0.0:Δr₊:Cell_R⁺ₛ)
Rs⁻ₛ = collect(0.0:Δr₋:Cell_R⁻ₛ)

with a simple for loop:

Rs⁺ₛ = zeros(Float64,50,1)
Rs⁻ₛ = zeros(Float64,50,1)
for i in 1:50
     Rs⁺ₛ[i] = Cell_R⁺ₛ * (i-1)
     Rs⁻ₛ[i] = Cell_R⁻ₛ * (i-1)
end

And now it works!
But I still don't know why collect was causing such a catastrophic issue.

[ChrisRackauckas]

Honed it down a bit:

using DifferentialEquations
using DiffEqFlux, Optim
using Flux

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

    """
    #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.
    """
    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.
    """
    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.
    """
    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)
    @show typeof(Δr₊),typeof(Cell_R⁺ₛ)
    Rs⁺ₛ = collect(0.0:Δr₊:Cell_R⁺ₛ)
    Rs⁻ₛ = collect(0.0:Δr₋:Cell_R⁻ₛ)
    dCₜ .= 0
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)

function loss(p)
    tmp_prob = remake(SPM_model,p=p)
    tmp_sol = solve(tmp_prob,alg_hints = [:stiff])
    sum(tmp_sol)
end

using ForwardDiff
ForwardDiff.gradient(loss, p)

This led me to the MWE:

using ForwardDiff
function f(p)
    sum(collect(0.0:p[1]:p[2]))
end
ForwardDiff.gradient(f, [0.2,25.0])

[ChrisRackauckas]

Opened an upstream issue: JuliaDiff/ForwardDiff.jl#548

[ChrisRackauckas]

Fixed upstream.