I am trying to fit a differential equation (using DiffEqParamEstim) with these properties:

- Lumped parameter linear thermal model with 3 states
- 8 interpolated external inputs (sampled at 50 Hz) using Interpolations.jl objects (shown below with
`_itp`

suffix) - 8 parameters to fit
- Long simulation time (data length - 600s) - runs in 0.01 secs

I have pasted a snippet of the code below for the differential equation. The variables not mentioned in the function are all constant floats enclosed in a wrapper function. (apologies for the incomplete code, I tried to make a standalone example but couldnât replicate the slow behavior)

```
function thermal_model!(du, u, p, t)
# scaling parameters
p_friction, p_forcedconv, p_tread2road, p_deflection,
p_natconv, p_carcass2air, p_tread2carcass, p_air2ambient = p
fxtire = fxtire_itp(t)
fytire = fytire_itp(t)
fztire = fztire_itp(t)
vx = vx_itp(t)
alpha = alpha_itp(t)
kappa = kappa_itp(t)
r_loaded = r_loaded_itp(t)
h_splitter = h_splitter_itp(t)
# arc length of tread area
theta_1 = acos(min(r_loaded - h_splitter, r_unloaded) / r_unloaded)
theta_2 = acos(min(r_loaded, r_unloaded) / r_unloaded)
area_tread_forced_air = r_unloaded * (theta_1 - theta_2) * tire_width
area_tread_contact = tire_width * 2 * sqrt(max(r_unloaded^2 - r_loaded^2, 0))
q_friction = p_friction * vx * (abs(fytire * tan(alpha)) + abs(fxtire * kappa))
q_tread2ambient_forcedconv = p_forcedconv * h_forcedconv * area_tread_forced_air * (t_tread - t_ambient) * vx^0.805
q_tread2ambient_natconv = p_natconv * h_natconv * (area_tread - area_tread_contact) * (t_tread - t_ambient)
q_tread2carcass = p_tread2carcass * h_tread2carcass * area_tread * (t_tread - t_carcass)
q_carcass2air = p_carcass2air * h_carcass2air * area_tread * (t_carcass - t_air)
q_carcass2ambient_natconv = p_natconv * h_natconv * area_sidewall * (t_carcass - t_ambient)
q_tread2road = p_tread2road * h_tread2road * area_tread_contact * (t_tread - t_track)
q_deflection = p_deflection * h_deflection * vx * abs(fztire)
q_air2ambient = p_air2ambient * h_natconv * area_rim * (t_air - t_ambient)
du[1] = der_t_tread = (q_friction - q_tread2carcass - q_tread2road - q_tread2ambient_forcedconv - q_tread2ambient_natconv)/(m_tread * cp_tread)
du[2] = der_t_carcass = (q_tread2carcass + q_deflection - q_carcass2air - q_carcass2ambient_natconv)/(m_carcass * cp_carcass)
du[3] = der_t_air = (q_carcass2air - q_air2ambient)/(m_air * cp_air)
end
```

I first use BlackBoxOptim with parameter bounds to first get the model to a reasonable parameter space and then use that as an initial guess for BFGS(). I find that BlackBoxOptim is able to get the fit close but BFGS() is not refining the fit even after running for a long time. On top of that, it is also extremely slow.

Loss Function:

```
cost_function = build_loss_objective(prob, Tsit5(), L2Loss(tsteps, transpose(data_ref)),maxiters=500000,verbose=true)
```

This is the code for BFGS() optimization which uses the result from BlackBoxOptim.

```
n_p = 8
l_b = repeat([0.0],n_p)
u_b = repeat([2.0e2],n_p)
result_bfgs = Optim.optimize(cost_function,l_b,u_b,p_bb, Fminbox(BFGS()))
```

My guess is it is something to with the gradient calculation but I donât know how I can make it faster. I am using the default Tsit5() solver with adaptive timestep. I tried these approaches and failed:

- swapped to DataInterpolations.jl and it was even slower.
- tried fitting only 100s but there is drift in the states (its still relatively slower than bboptimize)

I am new to Julia and I would appreciate any pointers on what might be the issue here.