MaximeVH / EquivalentCircuits.jl

A julia package to either fit the parameters of a specified equivalent electrical circuit to electrochemical impedance data, or to suggest a plausible circuit configuration for a given set of measurements (either through a comparison of circuits from the literature, or through an evolutionary algorithm approach).
MIT License
22 stars 6 forks source link

Use Optimization.OptimizationProblem() to find optimal fit of parameters #19

Closed StefanPofahl closed 1 year ago

StefanPofahl commented 1 year ago

Hi Maxime,

could you help me to understand, what needs to be done to use the output of your function EquivalentCircuits.objectivefunction() in the optimization package Optimization.jl? Please have a look on my code: EC_optimisation_via_evolution_strategy.jl

Regards,

Stefan

MaximeVH commented 1 year ago

The Optimization.jl package requires the function to be optimized to have two inputs, but the objectivefunction() outputs a function that only takes one input. You can solve this as follows:

  1. Implement a version of objectivefunction() that outputs a function with two inputs (one of which will be unused):

    function objectivefunction_2input(circuitfunc,measurements,frequencies,weights=nothing) 
    if isnothing(weights) ; weights = ones(length(frequencies)) ; end
    function objective(x,p)
        model_output = [circuitfunc(x,fr) for fr in frequencies] 
        return  mean(weights .* (abs.(measurements - model_output).^2)./(abs.(measurements).^2 .+ abs.(model_output).^2))
    end
    return objective
    end
  2. Adjust the function that does the parameter optimisation to the syntax of Optimization.jl:

    
    using Optimization, OptimizationCMAEvolutionStrategy, OptimizationEvolutionary

function parameopt(circuitstring::String,measurements,frequencies,method) elements = foldl(replace,["["=>"","]"=>"","-"=>"",","=>""],init = denumber_circuit(circuitstring)) initial_parameters = EquivalentCircuits.flatten(EquivalentCircuits.karva_parameters(elements)); circfunc = circuitfunction(circuitstring) objective = objectivefunction_2input(circfunc,measurements,frequencies) lower = zeros(length(initial_parameters)) upper = get_parameter_upper_bound(circuitstring)

prob = OptimizationProblem(objective, init_params,SciMLBase.NullParameters(),lb = lower , ub = upper)
parameters = solve(prob, method)

return parameters

end

After these adjustments, it should work, for example:
```julialang
method_ = Evolutionary.CMAES(μ = 40, λ = 100)
parameopt("R1-L2-[P3,R4]-[P5,R6]-[P7,R8]",Z_data,frequ_data,method_)

Note that I removed the second optimisation fine-tuning step for clarity. You can copy-paste it from the original parameteroptimisation function if you need it.

Regards,

Maxime

StefanPofahl commented 1 year ago

Hi Maxime,

thanks for the quick response :-) Another tip came from Fredrik Bagge Carlson:

obj = (x, p) -> objective(x)

https://discourse.julialang.org/t/equivalent-circuit-optimization/93907/8

P.S.: Unfortunately, the package Optimization.jl seem to be quite buggy. I have not manged to use the variant ".ES": https://docs.sciml.ai/Optimization/stable/optimization_packages/evolutionary/ https://wildart.github.io/Evolutionary.jl/stable/es/

Following call crashes on my machine:

sol = solve(prob, Evolutionary.ES(μ = 40, λ = 100))

or:

sol = solve(prob, Evolutionary.ES(μ = 40, ρ = 3, λ = 100)) # three parents have one offspring