How to calculate quality of fit result?

[StefanPofahl]

I spent some time to figure out, how to compute a quality criteria by means of your functions, and I have the impression I do not understand the meaning of the function: circuitfitness(). Below my sample code, which includes my own approach to calculate a quality criteria. Does my code make sense? What would be the correct way, if one would like to use your functions?

using EquivalentCircuits, PlotlyJS, Printf, RobustModels

# --- sample impedance data:
impedance_values_ = ComplexF64[5919.90084073586 - 15.794826681804063im, 5919.575521325405 - 32.677443741883025im, 5918.183674897797 - 67.57666460870544im, 
5912.242152808868 - 139.49441081175667im, 5887.119965779756 - 285.73699600024963im, 5785.038233646888 - 566.878749499386im, 5428.935296370544 - 997.1881947423717im, 
4640.2144606930815 - 1257.8277219098052im, 3871.8361085209845 - 978.9656717819273im, 3537.682636142598 - 564.9627167404748im, 3442.9419240480606 - 315.3996363823805im, 
3418.140460121871 - 219.68986957046025im, 3405.513645508888 - 242.57272592660013im, 3373.904450003017 - 396.0671717029891im, 3249.673719526995 - 742.0301829777005im, 
2808.423185495856 - 1305.924162464249im, 1779.4087896271944 - 1698.9660879948128im, 701.9588433822435 - 1361.4674351816855im, 208.28978681589047 - 777.6453690080142im, 
65.93273498232111 - 392.50667235657im]

# --- corresponding frequencies:
frequency_values = [0.1, 0.20691380811147891, 0.42813323987193935, 0.8858667904100828, 1.8329807108324356, 3.792690190732248, 7.847599703514613, 16.237767391887218, 
33.59818286283781, 69.51927961775601, 143.84498882876616, 297.63514416313194, 615.8482110660267, 1274.2749857031336, 2636.650898730358, 5455.5947811685255, 11288.378916846883, 
23357.214690901215, 48329.30238571752, 100000.0]

# --- generate frequency vector with n_elements with the same range as given in the measurement:
n_elements = 100
frequ_vec = exp10.(LinRange(log10(frequency_values[1]), log10(frequency_values[end]), n_elements))

# --- examples of suitable equivalent circuits
nr_best_circuits = 4
if nr_best_circuits == 3
    circuit_model_preset = "[C1,[C2-[R3,C4],R5]]"
    circuit_params_preset   = (C1 = 2.322248710116646e-9, C2 = 7.146778669252158e-7, R3 = 8015.389370331851, C4 = 1.6325663887245989e-9, R5 = 5918.9481528813185)
elseif  nr_best_circuits == 2
    circuit_model_preset = "[C1,R2-[R3,C4]]"
    circuit_params_preset   = (C1 = 0.025036871360843482, R2 = 396.73873944116787, R3 = 2178.061127814435, C4 = 1.1589755057609664e-5)
elseif  nr_best_circuits == 3
    circuit_model_preset = "R1-[C2,R3-[C4,R5]]"
    circuit_params_preset   = (R1 = 20.012355936798915, C2 = 4.000335253194046e-9, R3 = 3400.1181419153604, C4 = 4.0010158529883644e-6, R5 = 2499.9496058123705)
    # circuit_par_preset   = (R1 = 18.133936476549355, C2 = 3.967856543272228e-9, R3 = 3401.2083814207344, C4 = 3.9972681328104986e-6, R5 = 2500.4785300668846)        
elseif  nr_best_circuits == 4
    circuit_model_preset = "[[C1,R2],[C3,R4]-R5]"
    circuit_params_preset   = (C1 = 3.95177182904296e-9, R2 = 12300.52241056916, C3 = 2.075152378369779e-6, R4 = 6687.854083227787, R5 = 4720.304069423031)
else
    error(string("Choise nr_best_circuits = ", nr_best_circuits,"does not exist!"))
end

# --- simulate impedance based on suitable preset of a circuit model and its corresponding parameter-set: 
circfunc_preset = EquivalentCircuits.circuitfunction(circuit_model_preset)
impedance_preset = EquivalentCircuits.simulateimpedance_noiseless(circfunc_preset, circuit_params_preset, frequ_vec)

# **************************************************************************************************************************
# --- function "circuitevolution()" to find a suitable equivalent circuit model and its parameters:
# **************************************************************************************************************************
# # Arguments
# - `generations::Integer=10`: the maximum number of iterations of the evolutionary algorithm.
# - `population_size::Integer=30`: the number of individuals in the population during each iteration.
# - `terminals::String="RCLP"`: the circuit components that are to be included in the circuit identification.
# - `head::Integer=8`: a hyperparameter than controls the maximum considered complexity of the circuits.
# - `cutoff::Float64=0.8`: a hyperparameter that controls the circuit complexity by removing redundant components.
# - `initial_population::Array{Circuit,1}=nothing`:the option to provide an initial population of circuits
#   (obtained by using the loadpopulation function) with which the algorithm starts.
# -------------------------------------------------------------------------------------------------------------------------
terminals_      = "RC"
head_           = 6
equiv_circ_evo  = EquivalentCircuits.circuitevolution(impedance_values_, frequency_values, terminals=terminals_, generations=100, population_size=20, head=head_)
# --- simulate Impedance:
circfunc_evo    = EquivalentCircuits.circuitfunction(equiv_circ_evo.Circuit)
impedance_evo   = EquivalentCircuits.simulateimpedance_noiseless(circfunc_evo, equiv_circ_evo.Parameters, frequ_vec)
# --- Calc quality as the mean of the distances between measured and simulated impedance:
impedance_evo_data_pts = EquivalentCircuits.simulateimpedance_noiseless(circfunc_evo, equiv_circ_evo.Parameters, frequency_values)
Q_              = RobustModels.mean(abs.(impedance_values_ - impedance_evo_data_pts))
println("Q: ", Q_, ", Circuit Model: ", equiv_circ_evo.Circuit)
println(equiv_circ_evo.Parameters)
# --- calc fitness:
karva_str       = EquivalentCircuits.generatekarva(head_, terminals_)
karva_parameters = EquivalentCircuits.karva_parameters(karva_str)
circuit_struct  = EquivalentCircuits.Circuit(karva_str, karva_parameters, nothing)
circuit_parameters, circuit_struct.fitness, param_inds = EquivalentCircuits.circuitfitness(circuit_struct, impedance_values_, frequency_values)
println(string("fitness: ", circuit_struct.fitness))

# --- plot measured against simulated impedance:
function plot_nyquist_comp_preset_evo()
    s_pts_info = Vector{String}(undef, 0)
    for i_ndx in eachindex(frequency_values)
        push!(s_pts_info, @sprintf("#:%i, f=%.2fHz", i_ndx, frequency_values[i_ndx]))
    end
    # ---
    trace_impedance = PlotlyJS.scatter(; x= real(impedance_values_), y=  imag(impedance_values_), name = "measured impedance", text = s_pts_info, mode = "markers")
    trace_simulated_preset  = PlotlyJS.scatter(; x= real(impedance_preset), y= imag(impedance_preset), name = string("preset: ", circuit_model_preset))
    trace_simulated_auto    = PlotlyJS.scatter(; x= real(impedance_evo),    y= imag(impedance_evo),    name = string("evolution: ",   equiv_circ_evo.Circuit))
    # ---
    plt_layout = PlotlyJS.Layout(;
        title_text          = "Sample <-> Preset <-> Evolution",
        xaxis_title_text    = "z<sub>Real</sub>",
        xaxis_dtick         = 1000,
        yaxis_title_text    = "z<sub>Imag</sub>",
        yaxis_dtick         = 1000,
        # --- Fixed Ratio Axes Configuration
        yaxis_scaleanchor   = "x",
        yaxis_scaleratio    = 1,
    )
    return PlotlyJS.Plot([trace_impedance, trace_simulated_preset, trace_simulated_auto], plt_layout)
end

plot_nyquist_comp_preset_evo()

[MaximeVH]

In your code example, you used the mean of the absolute differences between the impedance measurements and the impedance values simulated from the circuit. There are many ways (objective functions) reported in the literature to evaluate and optimise the quality of the fit. I have compiled a table of them and conducted an evaluation of which one generally works best.

It turns out that an objective function based on modulus weighting (F6 in the table) typically works well when tested over a range of circuits. An appropriate weighting should be included in the objective function, as otherwise, the terms with the largest magnitude would dominate, causing information to be lost at frequencies where the impedance magnitudes are small. F6 is also the objective function used during parameter optimisation in the package.

[MaximeVH]

The circuitfitness() function uses an objective function and an optimisation algorithm to optimise the circuit parameters. The value of the objective function with the optimised parameters measures how well a given circuit fits the impedance measurements.

[StefanPofahl]

Hi Maxime,

good overview! Could you help me with the notation used above? What are \Theta, index i and index m? Does Eq. 7 already include the weighing option?

Stefan

[MaximeVH]

Hi Stefan

Sure, in the displayed objective functions $\Theta$ is the parameterset, $Z{t,n}(\Theta)$ represents the parameter-dependent model estimation of the impedance, while $Z{m,n}$ represents the measured value of the $n$ th measurement (the index $m$ stands for measured). Single primes are the real parts of the complex impedance, whereas double primes are the complex parts. F7 does another kind of weighing and F12 is an example where no weighing is used. Finally, $\delta$ represents the euclidean distance between two impedance spectra. There seems to be a typo in F13, this should have been $$F13(\Theta) = \sum_{n=1}^{N} \delta(f_n)(\Theta)|fn - f{n-1}|,$$ where $f_{0}=0$, so there is no index $i$.

[StefanPofahl]

Thanks :-)