Nelder-Mead or other optimization

Specific Domains Optimization (Mathematical)
, , ,
1

I am trying to use Optim to find parameters of a vector valued function using the Nelder-Mead algorithm, but open to considering other algorithms or packages. The problem is to find parameters of a function by minimizing the sum of squares. I am getting the error message below.
Any ideas on what is wrong with the code below or any other algorithms/packages that will be more suitable?

using Optim
function curvefunc(T, beta0, beta1, beta2, beta3, lambda0, lambda1)
    xp1 = (1 - exp.(-T ./ lambda0)) ./ (T ./ lambda0)
    xp2 = alpha1 - exp.(-T ./ lambda0)
    xp3 = (1 - exp.(-T ./ lambda1)) ./ (T ./ lambda1) - exp.(-T ./ lambda1)
    return beta0 + beta1*xp1 + beta2*xp2 + beta3*xp3
end

function curvfit(T, Data)
    result = Optim.optimize(sum((curvefunc(T, Data, beta0, beta1, beta2, beta3, lambda0, lambda1)-Data).^2),zeros(6), NelderMead())
    
end

T = [0.5, 1.0, 2.0, 3.0, 4.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0]
Data = [0.011, 0.013, 0.016, 0.019, 0.021, 0.026, 0.03, 0.035, 0.037, 0.038, 0.04]

This gives the following error message:

curvfit(T, Data)
ERROR: MethodError: no method matching curvefunc(::Vector{Float64}, ::Vector{Float64}, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64)

Closest candidates are:
  curvefunc(::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any)
2

you need to pass Optim.optimize a function which takes in a vector and returns the objective value. But even after fixing that, I couldn’t run your code because I don’t know what alpha1 is.

Here’s an example that should point you in the right direction:

using Optim

function curvefunc(t, beta0, beta1, beta2, beta3, lambda0, lambda1)
    alpha1 = 0.0 # TODO; what is alpha1?
    xp1 = (1 - exp(-t / lambda0)) / (t / lambda0)
    xp2 = alpha1 - exp(-t / lambda0)
    xp3 = (1 - exp(-t / lambda1)) / (t / lambda1) - exp(-t / lambda1)
    return beta0 + beta1*xp1 + beta2*xp2 + beta3*xp3
end

function curvfit(T, data)
    function objective_function(args)
        return sum((curvefunc(t, args...) - d)^2 for (t, d) in zip(T, data))
    end
    return Optim.optimize(objective_function, zeros(6), NelderMead())
end

T = [0.5, 1.0, 2.0, 3.0, 4.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0]
Data = [0.011, 0.013, 0.016, 0.019, 0.021, 0.026, 0.03, 0.035, 0.037, 0.038, 0.04]
result = curvfit(T, Data)
Optim.minimizer(result)
3

Thank you!
I mixed up the code I wrote for JuMP with that for Optim. alpha1 was supposed to be xp1. Using JuMP I get very different set of results, but I guess this is because Nelder-Mead returns local minimum. I have shown my code and results.

using JuMP, Optim, AmplNLWriter, Couenne_jll, DataFrames

function deterministic_global(Data,T)
    model = Model(() -> AmplNLWriter.Optimizer(Couenne_jll.amplexe))

    @variables(model, begin
    β0
    β1
    β2
    β3
    λ0
    λ1 
end)

    
    @NLexpression(model,xp1[t=1:length(T)],(1 - exp(-T[t] / λ0)) / (T[t] / λ0))
    @NLexpression(model,xp2[t=1:length(T)],xp1[t] - exp(-T[t] / λ0))
    @NLexpression(model,xp3[t=1:length(T)],(1 - exp(-T[t] / λ1)) / (T[t] / λ1) - exp(-T[t] / λ1))
    @NLobjective(model,Min, sum((β0 + β1*xp1[t] + β2*xp2[t] + β3*xp3[t]-Data[t])^2 for t in 1:length(T)))
    optimize!(model)

    beta0 = value.(β0)
    beta1 = value.(β1)
    beta2 = value.(β2)
    beta3 = value.(β3)
    
    lambda0 = value.(λ0)
    lambda1 = value.(λ1)

    return beta0, beta1, beta2, beta3, lambda0, lambda1
end

function curvefunc(t, beta0, beta1, beta2, beta3, lambda0, lambda1)
    xp1 = (1 - exp(-t / lambda0)) / (t / lambda0)
    xp2 = xp1  - exp(-t / lambda0)
    xp3 = (1 - exp(-t / lambda1)) / (t / lambda1) - exp(-t / lambda1)
    return beta0 + beta1*xp1 + beta2*xp2 + beta3*xp3
end

function curvfit(T, data)
    function objective_function(args)
        return sum((curvefunc(t, args...) - d)^2 for (t, d) in zip(T, data))
    end
    return Optim.optimize(objective_function, zeros(6), NelderMead())
end

T = [0.5, 1.0, 2.0, 3.0, 4.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0]
Data = [0.011, 0.013, 0.016, 0.019, 0.021, 0.026, 0.03, 0.035, 0.037, 0.038, 0.04]

result = curvfit(T, Data)
Optim.minimizer(result)

beta0, beta1, beta2, beta3, lambda0, lambda1 = deterministic_global(Data,T)
df = DataFrame("Parameters" =>     ["β0","β1", "β2", "β3", "λ0", "λ1"], "Global" => [beta0, beta1, beta2, beta3, lambda0, lambda1], "Nelder-Mead" => Optim.minimizer(result))


 Row │ Parameters  Global        Nelder-Mead 
     │ String      Float64       Float64
─────┼───────────────────────────────────────
   1 │ β0           0.0218042      0.0317701
   2 │ β1          -0.00582319    -0.0221575
   3 │ β2          -0.000650847   -0.0505105
   4 │ β3          -0.0221717     -0.11358
   5 │ λ0          -3.86381        0.0880731
   6 │ λ1          -8.96302        0.0664493
4

You will have to compare the objective values to see which is a better fit.

1 Like
5

I am finding it strange all methods except Nelder-Mead in Optim are failing. I have tried both derivative free methods (Simulated Annealing) and derivatives based methods. Code and output for simulated annealing is given below. Does Optim have a different way of setting out the problem for other methods or am I doing something incorrect?

using Optim
function curvefunc(t, beta0, beta1, beta2, beta3, lambda0, lambda1)
    xp1 = (1 - exp(-t / lambda0)) / (t / lambda0)
    xp2 = xp1  - exp(-t / lambda0)
    xp3 = (1 - exp(-t / lambda1)) / (t / lambda1) - exp(-t / lambda1)
    return beta0 + beta1*xp1 + beta2*xp2 + beta3*xp3
end

function curvfit(T, data)
    function objective_function(args)
        return sum((curvefunc(t, args...) - d)^2 for (t, d) in zip(T, data))
    end
    return Optim.optimize(objective_function, zeros(6), SimulatedAnnealing())
end

T = [0.5, 1.0, 2.0, 3.0, 4.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0]
Data = [0.011, 0.013, 0.016, 0.019, 0.021, 0.026, 0.03, 0.035, 0.037, 0.038, 0.04]
result = curvfit(T, Data)

curvfit (generic function with 1 method)

 * Status: failure (reached maximum number of iterations)

 * Candidate solution
    Final objective value:     8.562000e-03

 * Found with
    Algorithm:     Simulated Annealing

 * Convergence measures
    |x - x'|               = NaN ≰ 0.0e+00
    |x - x'|/|x'|          = NaN ≰ 0.0e+00
    |f(x) - f(x')|         = NaN ≰ 0.0e+00
    |f(x) - f(x')|/|f(x')| = NaN ≰ 0.0e+00
    |g(x)|                 = NaN ≰ 1.0e-08

 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    1000
    f(x) calls:    1001
6

Two likely things:

  1. you have no variable bounds, so when lambda0 = 0 the function is undefined.
  2. it seems like simulated annealing found a good solution, even if it reached the maximum number of iterations. In most cases, derivative free optimizers will terminate after a number of iterations, not based on convergence measures.