Problems with LsqFit on an ODE-based model

[Eben60]

I have a 1D experimental dataset and an ODE-based model. The model is able to describe experimental data reasonably well by manual fitting of the parameter:

The next step now is using LsqFit, at first trying to fit just one parameter:

julia> fit1 = curve_fit(model, tdata, ydata, [p0]); 
julia> fit1.converged
true
julia> plot(fit1.resid)

However it apparently doesn’t even try to fit:

julia> fit1.param[1] == p0
true

Now let’s change the initial value

julia> fit2 = curve_fit(model, tdata, ydata, [p0*1.5]);
julia> fit2.converged
true
julia> fit2.param[1] == p0*1.5
true # !!
julia> plot(fit2.resid)

The Jacobians are in both cases zero:

julia> all(fit1.jacobian .== 0)
true
julia> all(fit2.jacobian .== 0)
true

Now I’m completely at loss

[ChrisRackauckas]

Don’t use LsqFit.jl. We actively recommend against it because it’s generally not as stable as using Optimization.jl-based cost functions with a full optimizer.

[Eben60]

Tomorrow is a project discussion with my boss, and it would be much better if we have this way of getting our material properties in addition/comparison to those used till now. Thus I may follow the advice and switch to another package later on, but today just trying to get results, quick and dirty.

I don’t need high precision, nor even precision estimation. Providing my own numeric differentiation for the Jacobian worked: Fit now converges to the same value for any reasonable initial parameters. Still would be nice to understand, why it didn’t work without it.

function j_m(t,p)
    J = Array{Float64,2}(undef, length(t),length(p))
    δrel = 1e-4
    p0 = p[1]
    δpm = p0 * δrel/2
    J[:,1] = (model(t, [p0+δpm]) .- model(t, [p0-δpm])) ./ (2δpm)   
    return J
end