Slow decrease of loss and an error during fitting an ODE to data

New to Julia
, ,
4

Thank you for all your reply. I was able to debug the code and apparently the error was due to a dimension mismatch when I am calculating the loss. The loss should be loss = mean(abs2,sol[1,:].-T).

I am getting a good fit now. But during the initial iterations the loss is still decreasing slowly. In the final few iterations the loss jumps drastically and gives the solution. Is this the normal behaviour? This is how my loss function goes.

8.823401982574351
8.822957812054376
8.822513653748809
8.822069507105471
8.821625372889988
8.821181250556224
8.82073714046958
8.820293042239316
8.819848956197449
8.819404882133725
8.819404882133725
8.819533712044315
0.9350706276822296
0.18154377871977903
0.17983946840351075
0.17983946632392875
retcode: Success

I tried reducing the tolerance. My new solver settings is

Solver = Tsit5()
sol = solve(newprob,Solver,abstol = 1e-6,reltol = 1e-9,saveat = tdata)

Still the loss decreases slowly

8.83070912956643
8.83026475925042
8.82982040099848
8.829376054809883
8.828931720683078
8.828487398616993
8.828043088610425
8.827598790661725
8.827154504770306
8.826710230934976
8.826265969153756
8.825821719426086
8.825377481750598
8.824933256125755
8.824489042550749
8.82404484102433
8.823600651544869
8.823156474111531
8.82271230872314
8.82226815537842
8.821824014075686
8.821379884814268
8.820935767592546
8.820491662410072
8.820047569264565
8.820047569264565
8.820047580535817
0.9202558030635122
0.18140240144957576
0.17976875876525386
0.1797687582869255

The following is my MWE

using DifferentialEquations, Optimization, OptimizationPolyalgorithms, SciMLSensitivity
using ForwardDiff, Plots
using JLD2
using Statistics


# Define the ODE
function ode_model!(du,u,p,t,T∞)
    T = u[1]
    C₁ = p[1]
    du[1] = -(T-T∞)/C₁
end

# Generate data
tdata = 1650:1:8850
T∞ = 273.15
C₁_actual = 650.0
u0 = [315.15]
p_actual = [C₁_actual]
ode_actual(du,u,p,t) = ode_model!(du,u,p,t,T∞)
prob_actual = ODEProblem(ode_actual,u0,(tdata[1],tdata[end]),p_actual)
yactual = solve(prob_actual,saveat=tdata)

# Plotting actual values
plot(tdata,yactual[1,:],label="Actual data",xlabel="Time",ylabel="Temperature",title="Actual data")


# Parameter estimation
pguess = [1000.0]
ode_model1!(du,u,p,t) = ode_model!(du,u,p,t,T∞)
prob = ODEProblem(ode_model1!,u0,(tdata[1],tdata[end]),pguess)
initial_sol = solve(prob ,saveat = tdata)

# view guessed solution
plt = plot(tdata,initial_sol[1,:],label = ["T predicted"],color = :red)
plot!(tdata,yactual[1,:],label = ["T actual"],color = :black)

# Define cost function
function loss(newp)
    newprob = remake(prob,p = newp)
    sol = solve(newprob,saveat = tdata)
    loss = mean(abs2,sol[1,:].-yactual[1,:])
    return loss
end

function callback(state,l)
    display(l)
    newprob = remake(prob,p = state.u)
    sol = solve(newprob,saveat = tdata)
    plt = plot(sol[1,:],ylim = (240,360),label = ["T current prediction"],color = :red)
    plot!(plt,tdata,yactual[1,:],label=["T actual"],color = :black)
    display(plt)
    return false
end

adtype = AutoForwardDiff()
pguess = [1000.0]
optf = OptimizationFunction((x,_) -> loss(x),adtype)
optprob = OptimizationProblem(optf,pguess)
pfinal = solve(optprob,PolyOpt(),callback = callback,maxiters = 500)

I am getting a really good fit with this method. I just want to know more about this behaviour. So if you have any insights please do share. Thank you