forked from MaximeVH/EquivalentCircuits.jl
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathEC_optimisation_experiment_Maxime.jl
executable file
·252 lines (231 loc) · 16.7 KB
/
EC_optimisation_experiment_Maxime.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
using EquivalentCircuits, PlotlyJS, Printf, RobustModels
using Optim, BlackBoxOptim
# --- project pages of Optim, BlackBoxOptim: -------------------------------------------------------------------------------
# --- https://github.com/robertfeldt/BlackBoxOptim.jl
# --- https://julianlsolvers.github.io/Optim.jl/stable/
# --- ----------------------------------------------------------------------------------------------------------------------
# --- Selection of optimisation methods from BlackBoxOptim package (see below): --------------------------------------------
# --- remark: no. 14 may crash "simultaneous_perturbation_stochastic_approximation"
idx_methode_selection = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]
idx_methode_selection = [1, 4, 10, 12]
# idx_methode_selection = [1, 4, 10, 13, 14]
idx_methode_selection = [1, 4, 10] # fast
# --- select part of data points for equivalent circuit fitting:
idx_data = range(8, stop= 52)
# --- Measurements at low frequencies should be more important!
weight_tuning_low_frequ = [10.0, 10, 10, 10, 10, 10, 10, 10.0, 10, 10, 5, 5, 5, 5]
# --- measurement frequencies: --------------------------------------------------------------------------------------------
frequ_data = [0.0199553, 0.0251206, 0.0316296, 0.0398258, 0.0501337, 0.0631739, 0.0794492, 0.1001603, 0.1260081, 0.1588983,
0.2003205, 0.2520161, 0.316723, 0.400641, 0.5040323, 0.6334459, 0.7923428, 0.999041, 1.266892, 1.584686, 1.998082,
2.504006, 3.158693, 3.945707, 5.008013, 6.317385, 7.944915, 9.93114, 12.40079, 15.625, 19.86229, 24.93351, 31.25,
38.42213, 50.22321, 63.3446, 79.00281, 100.4464, 125.558, 158.3615, 198.6229, 252.4038, 315.5048, 397.9953, 505.5147,
627.7902, 796.875, 998.264, 1265.625, 1577.524, 1976.103, 2527.573]
# --- measured values impedance values (normed by max value): ------------------------------------------------------------
Z_data = ComplexF64[1.0 + 0.00033364003270441563im, 0.9959327788858627 - 0.004310226899902616im, 0.9942611041385668 - 0.005671195638022529im,
0.9896626966561299 - 0.00856668211620483im, 0.9848663962129538 - 0.011108044348854046im, 0.9839133853757104 - 0.011680892393098745im,
0.9803825583393658 - 0.014357655072569432im, 0.9780390890674556 - 0.016721955182452106im, 0.977945350296579 - 0.019502872051785466im,
0.9749352942095479 - 0.023195138082417213im, 0.973617743707785 - 0.027590444894622im, 0.9703941715315353 - 0.03213677528212766im,
0.9655562094124143 - 0.03752154689802784im, 0.9597964827130084 - 0.04379162912776074im, 0.9520005416017874 - 0.04925451638605793im,
0.9442670930044841 - 0.054795519286752116im, 0.9348307234029257 - 0.05974284330522906im, 0.9245767434109456 - 0.06472662129015795im,
0.9139634314639392 - 0.06913234352134902im, 0.9049228477838593 - 0.07387656686959375im, 0.8949136301380565 - 0.07923529993802828im,
0.8851127208719789 - 0.08469818719632546im, 0.8739630148471799 - 0.0921400040620134im, 0.8622300453591497 - 0.1000505147820834im,
0.8483046301745104 - 0.10995557823802357im, 0.8333168422532719 - 0.1212302692906579im, 0.8155116834962478 - 0.13342151721408374im,
0.7950141389312739 - 0.14564921910396153im, 0.7706472662129017 - 0.15724678814517012im, 0.7418694635538452 - 0.16734453685235623im,
0.708347437546544 - 0.1744582680199768im, 0.6755492831587884 - 0.17673924477796932im, 0.6430635913406207 - 0.17507798544965972im,
0.6145565895752072 - 0.17028168500648364im, 0.580763762674263 - 0.16053285283533744im, 0.5548293693984576 - 0.1501643032345084im,
0.5330351051696932 - 0.13955099128750204im, 0.5123240445155008 - 0.12800029163173163im, 0.4953156653109263 - 0.11759528806445062im,
0.4794894361612932 - 0.10735172349145677im, 0.4655952672336127 - 0.09793097701837802im, 0.45218020757930044 - 0.08854147680225807im,
0.44080136233680345 - 0.08018310306577858im, 0.4300318191050031 - 0.07178306765335403im, 0.4198247084984612 - 0.06336220140295694im,
0.41134134973414643 - 0.055706868448050506im, 0.40291527577425623 - 0.04717664029829761im, 0.39567135186930735 - 0.03873494320992798im,
0.3888857063997542 - 0.029465220312150108im, 0.38332908037058067 - 0.020200705123865372im, 0.37827760216224104 - 0.009806116975570635im,
0.3737000255177766 + 0.002913495258380507im]
# --- reference fit is computed via a commercial EIS-software: --------------------------------------------------------------
trace_ref_fit_SW = ComplexF64[0.9788651042462025 - 0.003427725046892952im, 0.9787042129564272 - 0.0041986283934134444im, 0.9785003527527543 - 0.005145223347737662im,
0.9782393462484965 - 0.006306332111776151im, 0.9779013997203362 - 0.007727852770328056im, 0.9774552153205819 - 0.009476179242838563im,
0.9768659237035731 - 0.011593226994444013im, 0.9760613725326219 - 0.01419779962088405im, 0.9749709741261627 - 0.017323133776195227im,
0.9734501264062743 - 0.02111374980931685im, 0.9713406718464188 - 0.025603088443209795im, 0.9684505314612848 - 0.030769668472705174im,
0.9645292196704338 - 0.036566648193514566im, 0.9591448671132816 - 0.043022435797319084im, 0.9523887936368954 - 0.04948223664885959im,
0.9442179647804041 - 0.05564232645454433im, 0.9350958290630905 - 0.06105262795583616im, 0.924987982855181 - 0.06586070781401777im,
0.9145015589339538 - 0.0701685599788112im, 0.9048243526644042 - 0.07415006287794139im, 0.8950106684802762 - 0.07881112049424521im,
0.8853693610889181 - 0.08439055008090053im, 0.8748621766817914 - 0.09155553424484979im, 0.8636873268852573 - 0.09990018322926661im,
0.8498691726172095 - 0.11037443576081168im, 0.8339703638271251 - 0.12180969637711742im, 0.8154128823748253 - 0.13378783963504204im,
0.7942784123186765 - 0.1454474185257693im, 0.7701271385576475 - 0.1562531797843441im, 0.7418846967259558 - 0.16570911888871886im,
0.7098260570782756 - 0.17258871607415407im, 0.6778320548678061 - 0.1755989237878614im, 0.6456734302156821 - 0.17484887477148447im,
0.6169363253992103 - 0.17095547213511444im, 0.5820365788624768 - 0.16196409353188074im, 0.5548571975533133 - 0.1514921236560258im,
0.5320564866437297 - 0.1401634302966748im, 0.5107034591263453 - 0.12731820019180443im, 0.49382330842857997 - 0.11564715504857762im,
0.4788316298662813 - 0.10435619961101911im, 0.4661622167087829 - 0.09449953735021457im, 0.4541975586774136 - 0.08541745893062654im,
0.44379749463126544 - 0.07803589286507319im, 0.43325665751158693 - 0.07102760719786393im, 0.4224787798478673 - 0.06387538861048644im,
0.41289733172108517 - 0.05686348014163662im, 0.4029795762886971 - 0.04810976004278965im, 0.3947009737016643 - 0.03873887736795084im,
0.3874833173425731 - 0.027874091160924586im, 0.38222507364956054 - 0.01707695883305848im, 0.3781417585287105 - 0.005417302648024877im,
0.37489234021412426 + 0.008160231785984899im]
# --- defined equavalent circuit model (notation according to Julia Package: "EquivalentCircuits")
circ_strg_ref = "R1-L2-[P3,R4]-[P5,R6]-[P7,R8]"
# --- Sample of available optimisation methods from BlackBoxOptim package:
bb_opt_methods = [:separable_nes, :xnes, :dxnes, :adaptive_de_rand_1_bin, :adaptive_de_rand_1_bin_radiuslimited,
:de_rand_1_bin, :de_rand_1_bin_radiuslimited, :de_rand_2_bin, :de_rand_2_bin_radiuslimited, :generating_set_search,
:probabilistic_descent, :resampling_memetic_search, :resampling_inheritance_memetic_search,
:simultaneous_perturbation_stochastic_approximation]
bb_opt_method_names = ["01#separable_nes", "02#xnes", "03#dxnes", "04#adaptive_de_rand_1_bin", "05#adaptive_de_rand_1_bin_radiuslimited",
"06#de_rand_1_bin", "07#de_rand_1_bin_radiuslimited", "08#de_rand_2_bin", "09#de_rand_2_bin_radiuslimited", "10#generating_set_search",
"11#probabilistic_descent", "12#resampling_memetic_search", "13#resampling_inheritance_memetic_search",
"14#simultaneous_perturbation_stochastic_approximation"]
bb_opt_methods = bb_opt_methods[idx_methode_selection]
bb_opt_method_names = bb_opt_method_names[idx_methode_selection]
# --- function "denumber_circuit()"
denumber_circuit(Circuit) = replace(Circuit, r"[0-9]" => "")
# --- fill /build weighting vector: ----------------------------------------------------------------------------------------
function weighting_vector(_Z_measured::Vector{ComplexF64}, _weight_tuning::Vector{Float64}=Vector{Float64}(undef, 0);
b_tail::Bool=true)
_n_measuremnt = length(_Z_measured); _n_weighing = length(_weight_tuning)
if isempty(_weight_tuning)
_weight_vec = ones(_n_measuremnt, )
elseif _n_weighing < _n_measuremnt
fill_vec = ones(_n_measuremnt - _n_weighing, )
if b_tail
_weight_vec = vcat(fill_vec, _weight_tuning)
else
_weight_vec = vcat(_weight_tuning, fill_vec)
end
elseif length(_weight_tuning) > _n_measuremnt
@warn("Weighing factor is larger then number of measured points! - Truncated weighing vector will be used!")
_weight_vec = _weight_tuning[1:_n_measuremnt]
end
return _weight_vec
end
# Function used to optimise the circuit parameters, using a method from BlackBoxOptim, followed by Nelder--Mead simplex fine-tuning.
function parameopt(_circuitstring::String, _Z_measured, _frequencies, _optim_method; weighting=nothing)
elements = foldl(replace,["["=>"","]"=>"","-"=>"",","=>""], init = denumber_circuit(_circuitstring))
initial_parameters = EquivalentCircuits.flatten(EquivalentCircuits.karva_parameters(elements));
circfunc = circuitfunction(_circuitstring)
objective = EquivalentCircuits.objectivefunction(circfunc, _Z_measured, _frequencies, weighting)
lower = zeros(length(initial_parameters))
upper = get_parameter_upper_bound(_circuitstring)
### First step ###
SR = Array{Tuple{Float64, Float64}, 1}(undef, length(initial_parameters))
for (e, (l, u)) in enumerate(zip(lower, upper))
SR[e] = (l, u)
end
res = bboptimize(objective; SearchRange = SR, Method = _optim_method, TraceMode = :silent); #MaxSteps=70000,
initial_parameters = best_candidate(res);
if _optim_method == :simultaneous_perturbation_stochastic_approximation
initial_parameters = clamp.(initial_parameters, 0.0, 1.0)
end
fitness_1 = best_fitness(res);
### Second step ###
inner_optimizer = NelderMead()
println("DBG: [min..max] of initial_parameters: [", minimum(initial_parameters), " .. ", maximum(initial_parameters), "], _optim_method: :", _optim_method)
results = optimize(objective, lower, upper, initial_parameters, Fminbox(inner_optimizer),
Optim.Options(time_limit = 50.0)); #20.0
fitness_2 = results.minimum
best = results.minimizer
parameters = fitness_2 < fitness_1 ? best : initial_parameters
# ---
return parameters
end
# Optional: The parameter bounds for this particular application can be made more specific than the bounds
# that are generally used in the package, which is application-agnostic.
function get_parameter_upper_bound(readablecircuit::String)
elements = foldl(replace,["["=>"","]"=>"","-"=>"",","=>""], init = denumber_circuit(readablecircuit))
ranges = Dict('R'=>400,'C'=>0.01,'L'=>5,'P'=>[400,1],'W'=>400,'+'=>0,'-'=>0)
return EquivalentCircuits.flatten([ranges[e] for e in elements])
end
# Optimise the circuit parameters and simulate the resulting impedance spectrum.
function optimise_and_simulate(_circuitstring, _Z_measured, _frequencies, _optim_method)
_parameters = parameopt(_circuitstring, _Z_measured, _frequencies, _optim_method; weighting = weight_vec)
_trace = simulateimpedance_noiseless(circuitfunction(_circuitstring), _parameters, _frequencies)
return _trace, _parameters
end
#Simulate impedance spectra for circuit parameters using the considered optimisation methods.
function generate_traces(_circuitstring, _ref_data, _Z_measured, _frequencies, _opt_methods, _method_names)
_names = Vector{String}(undef, length(_opt_methods) + 2)
_fited_par = Dict{Int, Union{Missing, Any}}()
_traces = Vector{Vector{ComplexF64}}(undef, length(_opt_methods) + 2)
_traces[1] = _Z_measured; _names[1] = "Z_measured"; _fited_par[1] = missing
_traces[2] = _ref_data; _names[2] = "Ref_SW_fit"; _fited_par[2] = missing
for _i in 3:length(_traces)
@info(string("--- generate_traces, #: ", _i - 2, ", next methode: ", _method_names[_i - 2], "\t-----------------------------"))
_traces[_i], _fited_par[_i] = optimise_and_simulate(_circuitstring, _Z_measured, _frequencies, _opt_methods[_i - 2])
_names[_i] = _method_names[_i - 2]
end
return _traces, _names, _fited_par
end
# The fitting errors calculated using the modulus weighted objective function,
# you can adjust the function to see other fitting quality metrics (e.g. removal of the denominator gives the MSE).
function trace_quality(_Z_measured::Vector{ComplexF64}, trace::Vector{ComplexF64};
weighting_fac::Vector{Float64}=Vector{Float64}(undef, 0))
_n_measuremnt = length(_Z_measured); _n_weighting = length(weighting_fac)
if isempty(weighting_fac)
weighting_fac = ones(_n_measuremnt, )
elseif _n_weighing != _n_measuremnt
@warn("Length of weighing factor non-equal to number of measured points! - Unity vector will be used!")
weighting_fac = ones(_n_measuremnt, )
end
return mean((abs.(weighting_fac .* (_Z_measured - trace)).^2)./(abs.(_Z_measured).^2 .+ abs.(trace).^2))
end
# --- Nyquist plots of arbitraty number of impedance spectra, compared to Z_measured and reference fit.
function plot_Nyquist(_traces::Vector{Vector{ComplexF64}}, names::Vector{String}, _frequ_data::Vector{Float64})
_dtick = 1.0e-1
s_pts_info = []
for i_ndx in eachindex(_frequ_data)
push!(s_pts_info, @sprintf("#:%i, f=%.2fHz", i_ndx, _frequ_data[i_ndx]))
end
# --- traces:
T = Array{GenericTrace{Dict{Symbol, Any}}}(undef,length(_traces))
for i in 1:length(_traces)
if i == 1
T[i] = PlotlyJS.scatter(; x= real(_traces[i]), y= imag(_traces[i]), name = names[i], text = s_pts_info, mode = "markers")
else
T[i] = PlotlyJS.scatter(; x= real(_traces[i]), y= imag(_traces[i]), name = names[i])
end
end
# --- layout:
plt_layout = PlotlyJS.Layout(;
title_text = "Comparison of Fit Results",
xaxis_title_text = "normed z<sub>Real</sub> / -",
xaxis_dtick = _dtick,
yaxis_title_text = "normed z<sub>Imag</sub> / -",
yaxis_dtick = _dtick,
yaxis_scaleanchor = "x",
yaxis_scaleratio = 1,
)
return PlotlyJS.Plot(T, plt_layout)
end
# --- Plot the fitting errors to compare the different optimisation methods. -----------------------------------------------
function Plot_fitting_errors(_method_names, fit_errors)
bar_plot = PlotlyJS.bar(x= _method_names, y= fit_errors, text= string.(fit_errors), textposition= "auto")
plt_layout = PlotlyJS.Layout(;
title_text = "Comparison of fitting methods",
yaxis_title_text = "Fitting error",
yaxis_scaleanchor = "x",
yaxis_scaleratio = 1,
)
return PlotlyJS.Plot(bar_plot,plt_layout)
end
# --- main part: -----------------------------------------------------------------------------------------------------------
n_measured = length(frequ_data)
if idx_data[end] > n_measured
error("Selection of measured points \"idx_data\" exceeds vector of measured points!")
end
# --- select data:
frequ_data_sel = frequ_data[idx_data]
Z_data_sel = Z_data[idx_data]
trace_ref_fit_SW_sel = trace_ref_fit_SW[idx_data]
# --- build weighting vector of the correct length:
weight_vec = weighting_vector(Z_data_sel, weight_tuning_low_frequ, b_tail= false)
# --- Evaluate the optimisation methods and make plots: --------------------------------------------------------------------
@time begin
traces_, names_, fitted_params = generate_traces(circ_strg_ref, trace_ref_fit_SW_sel, Z_data_sel, frequ_data_sel, bb_opt_methods, bb_opt_method_names)
end
# --- list qualities / deviation from measured values and sort fit results according to their quality: ---------------------
Fitting_err = [trace_quality(Z_data_sel, traces_[i]) for i in 1:length(traces_)]
idx_sorted = sortperm(Fitting_err[3:end]); idx_sorted = 2 .+ idx_sorted
Fitting_err_sorted = vcat(Fitting_err[1:2], Fitting_err[idx_sorted])
names_sorted = vcat(names_[1:2], names_[idx_sorted])
traces_sorted = vcat(traces_[1:2], traces_[idx_sorted])
# --- plot:
quality_plot = Plot_fitting_errors(names_sorted, Fitting_err_sorted)
fn_plt = raw"c:\tmp\plt\quality_plot.html"
PlotlyJS.savefig(quality_plot, fn_plt)
nyqs_plot = plot_Nyquist(traces_sorted[1:4], names_sorted[1:4], frequ_data_sel)
fn_plt = raw"c:\tmp\plt\nyquist_plot.html"
PlotlyJS.savefig(nyqs_plot, fn_plt)