Hello @ChrisRackauckas, full code below. The code needs input data (which I don’t own so can’t share in public), it won’t run without it.

```
## Functions and common constants for Thermal Tire Model
using Printf
using FileIO
using DataFrames
using Interpolations
using DifferentialEquations
using RecursiveArrayTools
## Build ODE Problem
function build_tire_thermal_ode(data,init_state,tspan,tstep,wheel_pos)
## convenience lambda functions to append wheel position
wheelvar = (var) -> @sprintf("%s%s",var,wheel_pos)
get_data = (var) -> data[!, wheelvar(var)]
lin_interp = (var) -> LinearInterpolation(data[!, "Time"], get_data(var))
get_init = (var) -> init_state[wheelvar(var)]
## Define Interpolations
# TODO: figure out how to do multiple 1D interpolations of matrix columns; this is painful
fxtire_itp = lin_interp("FxTire")
fytire_itp = lin_interp("FyTire")
fztire_itp = lin_interp("FzTire")
vx_itp = lin_interp("vContactPatchX")
alpha_itp = lin_interp("aTireSlip")
kappa_itp = lin_interp("rTireSlip")
r_loaded_itp = lin_interp("rTireLoaded")
h_splitter_itp = lin_interp("hSplitterDeflected")
t_carcass_itp = lin_interp("TTireCarcass")
t_tread_itp = lin_interp("TTireSurface")
t_air_itp = lin_interp("TTireAir")
## Nominal Tire Model parameters
# using parameters from Limebeer and Mavros papers and scaling them using the fitting process
# Environment
t_ambient = init_state["TAmbient"]
t_track = init_state["TTrack"]
tire_size = [305, 54, 15]
r_rim = 0.5 * tire_size[3] * 0.0254
r_unloaded = init_state[wheelvar("rTireUnloaded")]
tire_width = tire_size[1] / 1000
area_tread = 2 * pi * r_unloaded * tire_width
area_rim = 2 * pi * r_rim * tire_width
area_sidewall = pi * (r_unloaded^2 - r_rim^2)
# air
tire_air_volume = area_sidewall * tire_width
air_density = 2.5738
m_air = air_density * tire_air_volume
# tread
tread_depth = 0.00254
tread_density = 1.10e03
tread_volume = tread_depth * tread_density
m_tread = area_tread * tread_volume
# carcass
m_tire = 9.072
m_carcass = m_tire - m_tread
carcass_depth = tread_depth * 2
carcass_volume = carcass_depth * area_tread
# road
road_depth = 20e-3
# specific heat
cp_air = 0.718e3 # J/Kg/K
cp_carcass = 1.7652e3
cp_tread = 1.560e3
# heat transfer coeffs
h_tread2road = 7.5e-1 / road_depth
h_tread2carcass = 3.03e-1 / carcass_depth
h_natconv = 5.0
h_forcedconv = 7.0
h_carcass2air = 8.0e1
h_deflection = 2.79e-4*60.0
## Thermal Model - using nominal physical heat transfer coefficients
smoothabs = x -> sqrt(x^2 + 10^-10) # smooth version of abs
function tire_thermal_model!(du, u, p, t)
t_tread, t_carcass, t_air = u
# 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)
# crude approximation of tire contact area purely based on geometry
area_tread_forced_air = h_splitter * tire_width
area_tread_contact = tire_width * 2 * sqrt(r_unloaded^2 - r_loaded^2)
# heat transfer rates
q_friction = p_friction * vx * (smoothabs(fytire * tan(alpha)) + smoothabs(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 * smoothabs(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
## Build ODE Problem
# initial state
u0 = map(get_init,["TTireSurface","TTireCarcass","TTireAir"])
n_p = 8
p0 = repeat([0.1],n_p)
ode_prob = ODEProblem(tire_thermal_model!, u0, tspan, p0)
## Get Reference Data
tsteps = collect(tspan[1]:tstep:tspan[2])
data_ref = convert(Array, VectorOfArray([t_tread_itp(tsteps), t_carcass_itp(tsteps), t_air_itp(tsteps)]))
return (ode_prob,data_ref)
end
## ODE Solution Plotting
using Plots
plotly()
function solve_ode_prob(prob,p,tspan,tsteps)
prob_remade = remake(prob;p = p,tpsan = tspan)
return solve(prob_remade, Tsit5(), saveat = tsteps)
end
function evaluate_ode_prob(prob,p,tspan,tstep,data_ref,plot_title)
tsteps = collect(tspan[1]:tstep:tspan[2])
@time sol = solve_ode_prob(prob,p,tspan,tsteps)
plt = plot(sol, title = plot_title, label = ["Tread" "Carcass" "Air"], xlabel = "Time (s)", ylabel = "Temperature (K)",legend=:topleft)
# size=(1200,600)
plot!(tsteps, data_ref, label="", ls = :dot)
display(plt)
return sol
end
## Training - Build ODE Parameter Estimation Problem
using DifferentialEquations, DiffEqParamEstim
function build_ode_estimation_prob(tspan,tstep,data_train,init_train)
tsteps = collect(tspan[1]:tstep:tspan[2])
prob, data_ref = build_tire_thermal_ode(data_train,init_train,tspan,tstep,wheel_pos)
cost_function = build_loss_objective(prob, Tsit5(), L2Loss(tsteps, transpose(data_ref));
maxiters=500000,verbose=true)
return prob, data_ref,cost_function
end
## Load data
using DataFrames, Printf, FileIO
data_dir = joinpath(pwd(), "data", "processed")
wheel_pos = "RF"
wheelvar = (var) -> @sprintf("%s%s",var,wheel_pos)
filename_train = "Phoenix_R18_50Hz"
filepath_train = joinpath(data_dir,@sprintf("%s.jld2",filename_train))
filename_validate = "Phoenix_R17_50Hz"
filepath_validate = joinpath(data_dir,@sprintf("%s.jld2",filename_validate))
# read initial conditions and reference data
data_train,init_train = load(filepath_train, "df_data", "init_state")
data_validate,init_validate = load(filepath_validate, "df_data", "init_state")
## simulation constants
sampling_rate = 50.0 # Hz
tstep = 1/sampling_rate
t_end = 600.0 # end time
tspan = (0.0, t_end)
tsteps = collect(tspan[1]:tstep:tspan[2])
## Parameters
n_p = 8
param_names = [:p_friction,:p_forcedconv,:p_tread2road,:p_deflection, :p_natconv, :p_carcass2air, :p_tread2carcass, :p_air2ambient]
display_param = param_val -> display(Dict(param_names.=>param_val))
p0 = repeat([0.1], n_p)
## Build ODE Problem
prob, data_ref,cost_function = build_ode_estimation_prob((0.0,600.0),tstep,data_train,init_train)
data_ref_tr = transpose(data_ref)
## BlackBoxOptim
# [2.010205352273235, 3.6567658848983844, 86.61267821640898, 0.039428524429510406, 18.151228262926487, 0.0188954545517476, 16.84746373039449, 0.06875481085963224]
using BlackBoxOptim
res = bboptimize(cost_function; SearchRange = (0.0, 1.0e2),MaxSteps=100000, NumDimensions = 8)
p_bb = best_candidate(res)
display_param(p_bb)
_ = evaluate_ode_prob(prob,p_bb,(0.0,600.0),tstep,data_ref,"BBOptimize")
## BFGS
using Optim
prob, data_ref,cost_function = build_ode_estimation_prob((0.0,600.0),tstep,data_train,init_train)
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()))
print(result_bfgs)
p_bfgs = result_bfgs.minimizer
display_param(p_bfgs)
_ = evaluate_ode_prob(prob,p_bfgs,(0.0,600.0),tstep,data_ref,"BFGS")
## Validation
init_validate[wheelvar("TTireCarcass")] = init_train[wheelvar("TTireCarcass")]
init_validate[wheelvar("TTireAir")] = init_train[wheelvar("TTireAir")]
init_validate["TAmbient"] = 288.706
init_validate["TTrack"] = 301.483
prob_val, data_val = build_tire_thermal_ode(data_validate,init_validate,tspan,tstep,wheel_pos)
_ = evaluate_ode_prob(prob_val,p_bfgs,tspan,tstep,data_val,"Validation")
```

This is the profiling result (https://pastebin.com/7MWyzkrS) for BFGS optimization using the code below (I had to time limit the function otherwise it wouldn’t finish)

```
Profile.init(n = 10^8)
Profile.clear()
@profile Optim.optimize(cost_function,l_b,u_b,p_bb, Fminbox(BFGS()),Optim.Options(time_limit=300.0))
Profile.print(format=:flat)
```