Nonlinear curve fitting in Julia

Specific Domains Optimization (Mathematical)
,
1

Hi all,
My question is to do a nonlinear curve fitting task (bounded) for estimating 10 parameters from several give data points, presented as follows:

using Optim
# my data points
I = [100.254, 94.261, 88.224, 82.193, 76.215, 70.232, 64.205, 58.164, 52.103, 46.149, 40.147, 34.056, 28.049, 22.037, 16.035, 9.979, 3.932] # xdata
Es = [3.048, 3.124, 3.203, 3.279, 3.347, 3.415, 3.482, 3.544, 3.609, 3.674, 3.74, 3.811, 3.887, 3.971, 4.079, 4.207, 4.394] # ydata
# function to be fitted: p are the parameters array with 10 elements
function FCpol(I, p)
    # fitting parameters
    α_a, α_c, iloss_0, i0_a, i0_c = p[1], p[2], p[3], p[4], p[5]
    A1, Rion_0, R0, Bc_0, DO2_0 = p[6], p[7], p[8], p[9], p[10]
    Ncell = 5
    po2, ph2 = 0.21, 1
    Tfc, PO2 = 55 + 273.15, 130
    F, R, A0 = 96458, 8.314, 100
    Afc = A0*1 + A1*1 # ECSA
    i = I / Afc 
    Enerst = 1.229-8.5e-4(Tfc-298.15)+4.308e-5Tfc*(log(ph2)+0.5log(po2))
    Vact_a = (R*Tfc / (2α_a*F)) * log.((iloss_0*1 .+ I/Afc)/i0_a)
    Vact_c = (R*Tfc / (4α_c*F)) * log.((iloss_0 .+ i)/i0_c)
    Vohm = i .* (Rion_0 + R0)
    imax_c0 = 4F * (DO2_0) * PO2 
    Vcon = Bc_0 * log.(1 .-i/imax_c0)
    return Ncell * (Enerst.-Vact_a.-Vact_c.-Vohm.+Vcon)
end

function sqerror(p)
    err = 0.0
    for i in 1:length(I)
        pred_i = FCpol(I[i], p)
        err += (Es[i] - pred_i)^2
    end
    return err

end

lb = [0.5, 0.1, 1e-3, 0.0, 0.0, -1, 0.0, 0.0, 0.0, 0.1]
ub = [0.8, 0.5, 0.1, 0.01, 0.01, 1.0, 0.09, 0.09, 5.0, 0.9]
p0 = [0.6, 0.11, 0.01, 0.001, 0.001, 0.1, 0.01, 0.01, 1.0, 0.11]

res = optimize(sqerror, lb, ub, p0)
print(res)

But I run into errors:

Output exceeds the [size limit](command:workbench.action.openSettings?[). Open the full output data [in a text editor](command:workbench.action.openLargeOutput?c4aba7ce-b319-49f8-b3c3-35afb28528a8)

DomainError with -182400.38351824053: log will only return a complex result if called with a complex argument. Try log(Complex(x)). Stacktrace: [1] throw_complex_domainerror(f::Symbol, x::Float64) @ Base.Math ./math.jl:33 [2] _log(x::Float64, base::Val{:ℯ}, func::Symbol) @ Base.Math ./special/log.jl:304 [3] log @ ./special/log.jl:269 [inlined] [4] _broadcast_getindex_evalf @ ./broadcast.jl:670 [inlined] [5] _broadcast_getindex @ ./broadcast.jl:643 [inlined] [6] getindex @ ./broadcast.jl:597 [inlined] [7] copy @ ./broadcast.jl:875 [inlined] [8] materialize @ ./broadcast.jl:860 [inlined] [9] FCpol(I::Float64, p::Vector{Float64}) @ Main ~/Julia_projects/FCS_MTK.ipynb:20 [10] sqerror(p::Vector{Float64}) @ Main ~/Julia_projects/FCS_MTK.ipynb:30 [11] finite_difference_gradient!(df::Vector{Float64}, f::typeof(sqerror), x::Vector{Float64}, cache::FiniteDiff.GradientCache{Nothing, Nothing, Nothing, Vector{Float64}, Val{:central}(), Float64, Val{true}()}; relstep::Float64, absstep::Float64, dir::Bool)

...

@ ~/.julia/packages/Optim/6Lpjy/src/multivariate/solvers/constrained/fminbox.jl:269 [inlined] [27] optimize (repeats 3 times) @ ~/.julia/packages/Optim/6Lpjy/src/multivariate/solvers/constrained/fminbox.jl:265 [inlined] [28] top-level scope @ ~/Julia_projects/FCS_MTK.ipynb:41

In above solution I have refered to a solution in

Estimating parameters of a function for curve fitting

I have also tried the LsqFit.jl, but run into similar errors.
Thank you for checking my problem!

2

I think you should catch log domain error and/or return Inf when it occure or prevent it somehow.

3

Hi, thanks for the quick reply. It is weird, I have set the Bounds as in the program.
And also, I have tried Scipy curve_fit function, it works fine with the same data, without any errors.

4

Hi! Why are you using dot syntax with scalars? But this does’t matter much.

I think this is because there is “numerucal noise” with small numbers.

When you include in your function something like this:

if p[5] <= 0 p[5] = eps() end
if p[4] <= 0 p[4] = eps() end

It will works or change your bounds because i0_a, i0_c can’t be zero.

5

Thank you very much for the suggestion.
I tried to modify values in p0 of i0_a, i0_c into small numbers: 1e-5, but this time the programming never ends,
keeps running :joy:.