My lsqcurve fit is not working due YDATA and function value sizes are not being equal. (but they are?)

%I double checked the sizes and they are all set to 1,116. So I do not understand the problem.

This is my code:

Zreal = [101.8	104.7	105.7	106.8	106.7	108.5	107.7	108.9	108.1	108.6	109	108.8	109.1	107.8	109.5	109	109.9	109.9	108.4	110.7	110	112.4	110.7	112.1	112.9	114.1	115	116.2	118.8	121.3	125.8	127.3	130.8	122.8	138.3	143	146.9	151.9	158.3	162.5	169.1	177.9	189.7	205.3	228.2	259.1	304.5	366.1	453.5	568.1	730.6	933.1	1140	1368	1604	1746	1943	2126
];
Zim = [25.5800000000000	21.4300000000000	16.5700000000000	13.7500000000000	11.1300000000000	8.84500000000000	5.54700000000000	4.27600000000000	5.55100000000000	3.92100000000000	2.55200000000000	2.68300000000000	3.26400000000000	3.53800000000000	3.27400000000000	4.54300000000000	4.71500000000000	3.58100000000000	3.97600000000000	5.26300000000000	7.75600000000000	8.47000000000000	7.79900000000000	9.77800000000000	11.7600000000000	13.8800000000000	15.4500000000000	18.3700000000000	22.8200000000000	26.4400000000000	29.6200000000000	35.5500000000000	40.2700000000000	45.5900000000000	54.4900000000000	62.5100000000000	73.5100000000000	86.7100000000000	102.300000000000	125.200000000000	150.300000000000	181.400000000000	220.900000000000	270.200000000000	332	401.300000000000	486.300000000000	581.500000000000	689.500000000000	805.500000000000	915.900000000000	1008	1017	1010	999.300000000000	850.900000000000	818.700000000000	738.700000000000];
freq = [100100	79450	63140	50200	39890	31640	25170	20020	15890	12610	10080	8016	6328	5016	3984	3171	2528	1976	1578	1266	998.300000000000	796.900000000000	627.800000000000	505.500000000000	398	315.500000000000	252.400000000000	198.600000000000	158.400000000000	125.600000000000	100.400000000000	79	63.3400000000000	49.8700000000000	39.7200000000000	31.6700000000000	24.9300000000000	19.8600000000000	15.8400000000000	12.4000000000000	9.93100000000000	7.94500000000000	6.31700000000000	5.00800000000000	3.94600000000000	3.15900000000000	2.50400000000000	1.99800000000000	1.58500000000000	1.26700000000000	0.999000000000000	0.792300000000000	0.633400000000000	0.504000000000000	0.400600000000000	0.316700000000000	0.252000000000000	0.200300000000000];
parms = [2339	6.44000000000000e-05	68.0400000000000	0.00292000000000000	3.01500000000000	100.800000000000	0.000998600000000000	4.81200000000000	90.8500000000000	108.700000000000	0.000137600000000000	0.924800000000000];

concatenated_Z = [Zreal, Zim]';


size1 = size(my_fit_functioncpe(freq,parms))

size1 = 1×2

116 1

size2 = size(concatenated_Z)

size2 = 1×2

116 1

x = lsqcurvefit(@my_fit_functioncpe,[100 0],[freq,parms],concatenated_Z);

Error using lsqcurvefit
Function value and YDATA sizes are not equal.

function  concatz = my_fit_functioncpe(freq,parameters)
    Rsolc = parameters(1);
    Cm = parameters(2);
    Rq = parameters(3);
    Rn = parameters(4);
    Ln = parameters(5);
    RKch = parameters(6);
    Rm = parameters(7);
    Lm = parameters(8);
    Rl = parameters(9);
    Rctc = parameters(10);
    T = parameters(11);
    qc = parameters(12);

    w = 2 * pi * freq;
    Z = Rctc + 1./((1./Rsolc+T*(1j*w).^qc)) + 1./((1j*w*Cm)+ 1/Rl + 1/Rq + 1/RKch  + 1./(Rn + (1j*w*Ln)) + 1./(Rm + (1j*w*Lm)));%4c
    reZ = real(Z);
    imZ = imag(Z);
    concatz = cat(2,reZ,imZ)';
end

[Torsten]

Why do you still work with all these parameters to be fitted ?

It's like trying to fit a function of the form

y = (a+b)*x + c + d^2

You cannot get senseful results for a, b, c and d.

Say you get a = 2, b = 3, c = 4 and d = -1.

Is this better than

a = 3, b = 2, c = 1 and d = -2 ?

You see the problem ?

You can only fit two parameters, namely ab = a + b and cd = c + d^2:

y = ab*x + cd

Overfitting like above makes the Jacobian of the problem become singular.

Answer 1

[Star Strider]

There were several problems. First, the objective function arguments were reversed, and the data and calculations returned by it were definitely not the same sizes. This required a bit of tweaking, however it now works.

Try this —

Zreal = [101.8	104.7	105.7	106.8	106.7	108.5	107.7	108.9	108.1	108.6	109	108.8	109.1	107.8	109.5	109	109.9	109.9	108.4	110.7	110	112.4	110.7	112.1	112.9	114.1	115	116.2	118.8	121.3	125.8	127.3	130.8	122.8	138.3	143	146.9	151.9	158.3	162.5	169.1	177.9	189.7	205.3	228.2	259.1	304.5	366.1	453.5	568.1	730.6	933.1	1140	1368	1604	1746	1943	2126];
Zim = [25.5800000000000	21.4300000000000	16.5700000000000	13.7500000000000	11.1300000000000	8.84500000000000	5.54700000000000	4.27600000000000	5.55100000000000	3.92100000000000	2.55200000000000	2.68300000000000	3.26400000000000	3.53800000000000	3.27400000000000	4.54300000000000	4.71500000000000	3.58100000000000	3.97600000000000	5.26300000000000	7.75600000000000	8.47000000000000	7.79900000000000	9.77800000000000	11.7600000000000	13.8800000000000	15.4500000000000	18.3700000000000	22.8200000000000	26.4400000000000	29.6200000000000	35.5500000000000	40.2700000000000	45.5900000000000	54.4900000000000	62.5100000000000	73.5100000000000	86.7100000000000	102.300000000000	125.200000000000	150.300000000000	181.400000000000	220.900000000000	270.200000000000	332	401.300000000000	486.300000000000	581.500000000000	689.500000000000	805.500000000000	915.900000000000	1008	1017	1010	999.300000000000	850.900000000000	818.700000000000	738.700000000000];
freq = [100100	79450	63140	50200	39890	31640	25170	20020	15890	12610	10080	8016	6328	5016	3984	3171	2528	1976	1578	1266	998.300000000000	796.900000000000	627.800000000000	505.500000000000	398	315.500000000000	252.400000000000	198.600000000000	158.400000000000	125.600000000000	100.400000000000	79	63.3400000000000	49.8700000000000	39.7200000000000	31.6700000000000	24.9300000000000	19.8600000000000	15.8400000000000	12.4000000000000	9.93100000000000	7.94500000000000	6.31700000000000	5.00800000000000	3.94600000000000	3.15900000000000	2.50400000000000	1.99800000000000	1.58500000000000	1.26700000000000	0.999000000000000	0.792300000000000	0.633400000000000	0.504000000000000	0.400600000000000	0.316700000000000	0.252000000000000	0.200300000000000];
parms = [2339	6.44000000000000e-05	68.0400000000000	0.00292000000000000	3.01500000000000	100.800000000000	0.000998600000000000	4.81200000000000	90.8500000000000	108.700000000000	0.000137600000000000	0.924800000000000];

concatenated_Z = [Zreal; Zim]';


size1 = size(my_fit_functioncpe(parms,freq))

size1 = 1×2

58 2

size2 = size(concatenated_Z)

size2 = 1×2

58 2

x = lsqcurvefit(@my_fit_functioncpe,parms,freq,concatenated_Z)

Solver stopped prematurely.

lsqcurvefit stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 1.200000e+03.

x = 1×12

1.0e+03 * 1.4254 0.0000 0.0542 -0.0201 0.0001 0.0946 0.0933 0.0009 0.0832 0.1116 0.0000 0.0010

fittedZ = my_fit_functioncpe(x,freq);

figure
semilogx(freq, concatenated_Z, 'DisplayName','Z Data')
hold on
plot(freq, fittedZ, 'DisplayName','Fitted Z')
hold off
grid
xlabel('Frequency')
ylabel('Z')
legend('Location','best')

function  concatz = my_fit_functioncpe(parameters,freq)
    Rsolc = parameters(1);
    Cm = parameters(2);
    Rq = parameters(3);
    Rn = parameters(4);
    Ln = parameters(5);
    RKch = parameters(6);
    Rm = parameters(7);
    Lm = parameters(8);
    Rl = parameters(9);
    Rctc = parameters(10);
    T = parameters(11);
    qc = parameters(12);

    w = 2 * pi * freq;
    Z = Rctc + 1./((1./Rsolc+T*(1j*w).^qc)) + 1./((1j*w*Cm)+ 1/Rl + 1/Rq + 1/RKch  + 1./(Rn + (1j*w*Ln)) + 1./(Rm + (1j*w*Lm)));%4c
    reZ = real(Z);
    imZ = imag(Z);
    concatz = [reZ; imZ].';
    % concatz = cat(2,reZ,imZ)';
end

.

[Victor Carneiro da Cunha Martorelli]

Hi, Thank you so much for the reply! The graph seems a little off though. I am trying to optimise the parameters for the fit. However, the graph does not appear to be similar to the data. Also I need the minimum values of the parameters to be positive. So following your advice I got this, As the graph and function is supposed to fit the Z real and imaginary bits. Do you have any suggestions I could use to improve the fit? Thank you again so much.

Zreal = [101.8	104.7	105.7	106.8	106.7	108.5	107.7	108.9	108.1	108.6	109	108.8	109.1	107.8	109.5	109	109.9	109.9	108.4	110.7	110	112.4	110.7	112.1	112.9	114.1	115	116.2	118.8	121.3	125.8	127.3	130.8	122.8	138.3	143	146.9	151.9	158.3	162.5	169.1	177.9	189.7	205.3	228.2	259.1	304.5	366.1	453.5	568.1	730.6	933.1	1140	1368	1604	1746	1943	2126];
Zim = [25.5800000000000	21.4300000000000	16.5700000000000	13.7500000000000	11.1300000000000	8.84500000000000	5.54700000000000	4.27600000000000	5.55100000000000	3.92100000000000	2.55200000000000	2.68300000000000	3.26400000000000	3.53800000000000	3.27400000000000	4.54300000000000	4.71500000000000	3.58100000000000	3.97600000000000	5.26300000000000	7.75600000000000	8.47000000000000	7.79900000000000	9.77800000000000	11.7600000000000	13.8800000000000	15.4500000000000	18.3700000000000	22.8200000000000	26.4400000000000	29.6200000000000	35.5500000000000	40.2700000000000	45.5900000000000	54.4900000000000	62.5100000000000	73.5100000000000	86.7100000000000	102.300000000000	125.200000000000	150.300000000000	181.400000000000	220.900000000000	270.200000000000	332	401.300000000000	486.300000000000	581.500000000000	689.500000000000	805.500000000000	915.900000000000	1008	1017	1010	999.300000000000	850.900000000000	818.700000000000	738.700000000000];
freq = [100100	79450	63140	50200	39890	31640	25170	20020	15890	12610	10080	8016	6328	5016	3984	3171	2528	1976	1578	1266	998.300000000000	796.900000000000	627.800000000000	505.500000000000	398	315.500000000000	252.400000000000	198.600000000000	158.400000000000	125.600000000000	100.400000000000	79	63.3400000000000	49.8700000000000	39.7200000000000	31.6700000000000	24.9300000000000	19.8600000000000	15.8400000000000	12.4000000000000	9.93100000000000	7.94500000000000	6.31700000000000	5.00800000000000	3.94600000000000	3.15900000000000	2.50400000000000	1.99800000000000	1.58500000000000	1.26700000000000	0.999000000000000	0.792300000000000	0.633400000000000	0.504000000000000	0.400600000000000	0.316700000000000	0.252000000000000	0.200300000000000];
parms = [2339	6.44000000000000e-05	68.0400000000000	0.00292000000000000	3.01500000000000	100.800000000000	0.000998600000000000	4.81200000000000	90.8500000000000	108.700000000000	0.000137600000000000	0.924800000000000];
concatenated_Z = [Zreal; Zim]';

lb = zeros(1,12);
size1 = size(my_fit_functioncpe(parms,freq))

size1 = 1×2

58 2

size2 = size(concatenated_Z)

size2 = 1×2

58 2

x = lsqcurvefit(@my_fit_functioncpe,parms,freq,concatenated_Z,lb);

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

optfit = my_fit_functioncpe(x, freq);

figure()
scatter(Zreal,Zim)
hold on
plot(optfit(1:58),-optfit(59:116))
legend('Data','Optimized Fit')
hold off

function  concatz = my_fit_functioncpe(parameters,freq)
    Rsolc = parameters(1);
    Cm = parameters(2);
    Rq = parameters(3);
    Rn = parameters(4);
    Ln = parameters(5);
    RKch = parameters(6);
    Rm = parameters(7);
    Lm = parameters(8);
    Rl = parameters(9);
    Rctc = parameters(10);
    T = parameters(11);
    qc = parameters(12);

    w = 2 * pi * freq;
    Z = Rctc + 1./((1./Rsolc+T*(1j*w).^qc)) + 1./((1j*w*Cm)+ 1/Rl + 1/Rq + 1/RKch  + 1./(Rn + (1j*w*Ln)) + 1./(Rm + (1j*w*Lm)));%4c
    %CPEIMPEDANCE = 1./(T*(1j*w).^qc);
    %Z = 1./(1./(Rctc)+(1./(Rsolc+CPEIMPEDANCE))) + 1./((1j*w*Cm) +1./Rl+ 1./Rq + 1./RKch  + 1./(Rn + (1j*w*Ln)) + 1./(Rm + (1j*w*Lm)));
    reZ = real(Z);
    imZ = imag(Z);
    concatz = [reZ; imZ].';
end

[Star Strider]

My pleasure!

I usually rely on the ga (genetic algorithm) and occasionally other Global Optimization Toolbox functions to arrive at the best parameter estimates, so I let it do what it can to estimate the parameters. It is usually necessary to run it a number of times (sometimes in a loop) with different initial parameter estimates. You can then use the best estimates it returns (lowest fitness value) to determine which parameter set is the best. there is a 55 second limit here, so I need to constrain the number of generations so that it completes within that time. (Letting it run for a significantly longer time is generally preferable.) The 'InitialPopullationMatrix' is random, with its columns scaled by the magnitudes of your initial ‘parms’ vector.

Experiment with this, with increasing numbers of permitted generations. An options structure (that I usually use although will not work here) is:

opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms).*1E-3.*10.^round(log10(parms)), 'MaxGenerations',5E4, 'FunctionTolerance',1E-12, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);

You will simply need to change the options structure name in the ga call to use it.

Try this —

Zreal = [101.8	104.7	105.7	106.8	106.7	108.5	107.7	108.9	108.1	108.6	109	108.8	109.1	107.8	109.5	109	109.9	109.9	108.4	110.7	110	112.4	110.7	112.1	112.9	114.1	115	116.2	118.8	121.3	125.8	127.3	130.8	122.8	138.3	143	146.9	151.9	158.3	162.5	169.1	177.9	189.7	205.3	228.2	259.1	304.5	366.1	453.5	568.1	730.6	933.1	1140	1368	1604	1746	1943	2126];
Zim = [25.5800000000000	21.4300000000000	16.5700000000000	13.7500000000000	11.1300000000000	8.84500000000000	5.54700000000000	4.27600000000000	5.55100000000000	3.92100000000000	2.55200000000000	2.68300000000000	3.26400000000000	3.53800000000000	3.27400000000000	4.54300000000000	4.71500000000000	3.58100000000000	3.97600000000000	5.26300000000000	7.75600000000000	8.47000000000000	7.79900000000000	9.77800000000000	11.7600000000000	13.8800000000000	15.4500000000000	18.3700000000000	22.8200000000000	26.4400000000000	29.6200000000000	35.5500000000000	40.2700000000000	45.5900000000000	54.4900000000000	62.5100000000000	73.5100000000000	86.7100000000000	102.300000000000	125.200000000000	150.300000000000	181.400000000000	220.900000000000	270.200000000000	332	401.300000000000	486.300000000000	581.500000000000	689.500000000000	805.500000000000	915.900000000000	1008	1017	1010	999.300000000000	850.900000000000	818.700000000000	738.700000000000];
freq = [100100	79450	63140	50200	39890	31640	25170	20020	15890	12610	10080	8016	6328	5016	3984	3171	2528	1976	1578	1266	998.300000000000	796.900000000000	627.800000000000	505.500000000000	398	315.500000000000	252.400000000000	198.600000000000	158.400000000000	125.600000000000	100.400000000000	79	63.3400000000000	49.8700000000000	39.7200000000000	31.6700000000000	24.9300000000000	19.8600000000000	15.8400000000000	12.4000000000000	9.93100000000000	7.94500000000000	6.31700000000000	5.00800000000000	3.94600000000000	3.15900000000000	2.50400000000000	1.99800000000000	1.58500000000000	1.26700000000000	0.999000000000000	0.792300000000000	0.633400000000000	0.504000000000000	0.400600000000000	0.316700000000000	0.252000000000000	0.200300000000000];
parms = [2339	6.44000000000000e-05	68.0400000000000	0.00292000000000000	3.01500000000000	100.800000000000	0.000998600000000000	4.81200000000000	90.8500000000000	108.700000000000	0.000137600000000000	0.924800000000000];
concatenated_Z = [Zreal; Zim]';

lb = zeros(1,12);
size1 = size(my_fit_functioncpe(parms,freq))

size1 = 1×2

58 2

size2 = size(concatenated_Z)

size2 = 1×2

58 2

ftns = @(theta) norm(concatenated_Z - my_fit_functioncpe(theta,freq));
PopSz = 100;
Parms = 12;

optsAns = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms).*1E-3.*10.^round(log10(parms)), 'MaxGenerations',5E4, 'FunctionTolerance',1E-12);                                                                  % Options Structure For 'Answers' Problems
t0 = clock;
fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)

Start Time: 2023-12-21 03:59:36.8470

[x,fval,exitflag,output,population,scores] = ga(ftns, Parms, [],[],[],[],zeros(Parms,1),Inf(Parms,1),[],[],optsAns);

ga stopped because the average change in the fitness value is less than options.FunctionTolerance.

t1 = clock;
fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)

Stop Time: 2023-12-21 03:59:38.5553

GA_Time = etime(t1,t0);
DT_GA_Time = datetime([0 0 0 0 0 GA_Time], 'Format','HH:mm:ss.SSS');
fprintf('\nElapsed Time: %23.15E\t\t%s\n\n', GA_Time, DT_GA_Time)

Elapsed Time:   1.708307000000005E+00		00:00:01.708

fprintf('Fitness value at convergence = %.4f\nGenerations \t\t\t\t = %d\n\n',fval,output.generations)

Fitness value at convergence = 1875.9372
Generations 				 = 332

fprintf(1,'\tParameters:\n')

	Parameters:

for k1 = 1:length(x)
    fprintf(1, '\t\tParameter(%2d) = %8.5f\n', k1, x(k1))
end

		Parameter( 1) = 1982.25220
		Parameter( 2) =  0.00026
		Parameter( 3) = 323.04402
		Parameter( 4) =  0.00151
		Parameter( 5) = 21.31300
		Parameter( 6) = 876.69366
		Parameter( 7) =  0.00115
		Parameter( 8) = 103.74891
		Parameter( 9) = 830.53841
		Parameter(10) = 327.14997
		Parameter(11) =  0.00003
		Parameter(12) =  2.96400

% x = lsqcurvefit(@my_fit_functioncpe,parms,freq,concatenated_Z,lb);

optfit = my_fit_functioncpe(x, freq);

figure()
scatter(Zreal,Zim)
hold on
plot(optfit(1:58),-optfit(59:116))
legend('Data','Optimized Fit')
hold off

function  concatz = my_fit_functioncpe(parameters,freq)
    Rsolc = parameters(1);
    Cm = parameters(2);
    Rq = parameters(3);
    Rn = parameters(4);
    Ln = parameters(5);
    RKch = parameters(6);
    Rm = parameters(7);
    Lm = parameters(8);
    Rl = parameters(9);
    Rctc = parameters(10);
    T = parameters(11);
    qc = parameters(12);

    w = 2 * pi * freq;
    Z = Rctc + 1./((1./Rsolc+T*(1j*w).^qc)) + 1./((1j*w*Cm)+ 1/Rl + 1/Rq + 1/RKch  + 1./(Rn + (1j*w*Ln)) + 1./(Rm + (1j*w*Lm)));%4c
    %CPEIMPEDANCE = 1./(T*(1j*w).^qc);
    %Z = 1./(1./(Rctc)+(1./(Rsolc+CPEIMPEDANCE))) + 1./((1j*w*Cm) +1./Rl+ 1./Rq + 1./RKch  + 1./(Rn + (1j*w*Ln)) + 1./(Rm + (1j*w*Lm)));
    reZ = real(Z);
    imZ = imag(Z);
    concatz = [reZ; imZ].';
end

You will likely have to run it several times before it converges on good parameter estimates. I believe it can get there, although it is not getting there just yet. (In my experiments here, it continually converges on the negative of the desired result, however that is otherwise better than other results since it spans the appropriate ranges.)

in the loop, I would add a writematrix or writetable call that saves the necessary information (parameter estimates, fitness values, generations required) in each iteration, and do something else while you let it run for perhaps 100 or more iterations, then see what the best result is and the parameters that produced it. Then use these as the initial parameter estimates for lsqcurvefit to see if it can improve on them.

EDIT — (21 Dec 2023 at 05:06)

Thinking about it a few minutes later, since you are plotting:

plot(optfit(1:58),-optfit(59:116))

(that I did not pick up on until I came back to check it), that definitely explains the negative result in the plot. Your objective function (at least the one I am using with ga) probably needs to be negated as well in that respect. The ga function will likely converge on the correct parameters and the plot will match appropriately if you give it the same function you are using in the plot.

.

[Victor Carneiro da Cunha Martorelli]

Hi! Thank you again for the help. When you say to make a loop, do you mean to write a matrix file and then import it again after each iteration? I am not sure how that would fit the feed in the ga protocol. I understand what you mean but I am not sure if the parameters are being used in the ga sequence?

Sorry about this I really appreciate the help.

[Star Strider]

My pleasure!

I was thinking of something like this —

Zreal = [101.8 104.7 105.7 106.8 106.7 108.5 107.7 108.9 108.1 108.6 109 108.8 109.1 107.8 109.5 109 109.9 109.9 108.4 110.7 110 112.4 110.7 112.1 112.9 114.1 115 116.2 118.8 121.3 125.8 127.3 130.8 122.8 138.3 143 146.9 151.9 158.3 162.5 169.1 177.9 189.7 205.3 228.2 259.1 304.5 366.1 453.5 568.1 730.6 933.1 1140 1368 1604 1746 1943 2126];

Zim = [25.5800000000000 21.4300000000000 16.5700000000000 13.7500000000000 11.1300000000000 8.84500000000000 5.54700000000000 4.27600000000000 5.55100000000000 3.92100000000000 2.55200000000000 2.68300000000000 3.26400000000000 3.53800000000000 3.27400000000000 4.54300000000000 4.71500000000000 3.58100000000000 3.97600000000000 5.26300000000000 7.75600000000000 8.47000000000000 7.79900000000000 9.77800000000000 11.7600000000000 13.8800000000000 15.4500000000000 18.3700000000000 22.8200000000000 26.4400000000000 29.6200000000000 35.5500000000000 40.2700000000000 45.5900000000000 54.4900000000000 62.5100000000000 73.5100000000000 86.7100000000000 102.300000000000 125.200000000000 150.300000000000 181.400000000000 220.900000000000 270.200000000000 332 401.300000000000 486.300000000000 581.500000000000 689.500000000000 805.500000000000 915.900000000000 1008 1017 1010 999.300000000000 850.900000000000 818.700000000000 738.700000000000];

freq = [100100 79450 63140 50200 39890 31640 25170 20020 15890 12610 10080 8016 6328 5016 3984 3171 2528 1976 1578 1266 998.300000000000 796.900000000000 627.800000000000 505.500000000000 398 315.500000000000 252.400000000000 198.600000000000 158.400000000000 125.600000000000 100.400000000000 79 63.3400000000000 49.8700000000000 39.7200000000000 31.6700000000000 24.9300000000000 19.8600000000000 15.8400000000000 12.4000000000000 9.93100000000000 7.94500000000000 6.31700000000000 5.00800000000000 3.94600000000000 3.15900000000000 2.50400000000000 1.99800000000000 1.58500000000000 1.26700000000000 0.999000000000000 0.792300000000000 0.633400000000000 0.504000000000000 0.400600000000000 0.316700000000000 0.252000000000000 0.200300000000000];

parms = [2339 6.44000000000000e-05 68.0400000000000 0.00292000000000000 3.01500000000000 100.800000000000 0.000998600000000000 4.81200000000000 90.8500000000000 108.700000000000 0.000137600000000000 0.924800000000000];

concatenated_Z = [Zreal; Zim]';

lb = zeros(1,12);

size1 = size(my_fit_functioncpe(parms,freq))

size2 = size(concatenated_Z)

ftns = @(theta) norm(concatenated_Z - my_fit_functioncpe(theta,freq));

PopSz = 100;

Parms = 12;

optsAns = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms).*1E-3.*10.^round(log10(parms)), 'MaxGenerations',5E4, 'FunctionTolerance',1E-12); % Options Structure For 'Answers' Problems

writematrix(zeros(1,15),'Parameters.txt')

for k = 1:250

t0 = clock;

fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)

[x,fval,exitflag,output,population,scores] = ga(ftns, Parms, [],[],[],[],zeros(Parms,1),Inf(Parms,1),[],[],optsAns);

t1 = clock;

fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)

GA_Time = etime(t1,t0);

DT_GA_Time = datetime([0 0 0 0 0 GA_Time], 'Format','HH:mm:ss.SSS');

fprintf('\nElapsed Time: %23.15E\t\t%s\n\n', GA_Time, DT_GA_Time)

fprintf('Fitness value at convergence = %.4f\nGenerations \t\t\t\t = %d\n\n',fval,output.generations)

fprintf(1,'\tParameters:\n')

for k1 = 1:length(x)

fprintf(1, '\t\tParameter(%2d) = %8.5f\n', k1, x(k1))

end

addrow = [x(:).', fval, output.generations, GA_Time];

v = readmatrix('Parameters.txt');

writematrix([v; addrow],'Parameters.txt')

end

% x = lsqcurvefit(@my_fit_functioncpe,parms,freq,concatenated_Z,lb);

v = readmatrix('Parameters.txt')

vs = sortrows(v, 13, 'ascend')

optfit = my_fit_functioncpe(vs(2,1:12), freq);

figure()

scatter(Zreal,Zim)

hold on

plot(optfit(1:58),-optfit(59:116))

legend('Data','Optimized Fit')

hold off

function concatz = my_fit_functioncpe(parameters,freq)

Rsolc = parameters(1);

Cm = parameters(2);

Rq = parameters(3);

Rn = parameters(4);

Ln = parameters(5);

RKch = parameters(6);

Rm = parameters(7);

Lm = parameters(8);

Rl = parameters(9);

Rctc = parameters(10);

T = parameters(11);

qc = parameters(12);

w = 2 * pi * freq;

Z = Rctc + 1./((1./Rsolc+T*(1j*w).^qc)) + 1./((1j*w*Cm)+ 1/Rl + 1/Rq + 1/RKch + 1./(Rn + (1j*w*Ln)) + 1./(Rm + (1j*w*Lm)));%4c

%CPEIMPEDANCE = 1./(T*(1j*w).^qc);

%Z = 1./(1./(Rctc)+(1./(Rsolc+CPEIMPEDANCE))) + 1./((1j*w*Cm) +1./Rl+ 1./Rq + 1./RKch + 1./(Rn + (1j*w*Ln)) + 1./(Rm + (1j*w*Lm)));

reZ = real(Z);

imZ = imag(Z);

concatz = [reZ; imZ].';

end

The ‘Parameters.txt’ file then contains in each row:

[parameters, fitness_value, number_of_generations, elapsed_time]

If you want to, it would be relatively straightforward to choose the run with the best fitness by using sortrows on column 13:

v = readmatrix('Parameters.txt');

vs = sortrows(v, 13, 'ascend')

The matrix will then be sorted by fitness values, with the best fit in the first row. I chose to iterate 250 times in this example. Change that to whatever works best for you.

I tested this with two iterations and it works correctly.

.

[Victor Carneiro da Cunha Martorelli]

Thank you so much for the help this definitely works so much better!

[Star Strider]

As always, my pleasure!

[Alex Sha]

The result obtained below seems to be not bad:

Sum Squared Error (SSE): 7003.74699273406
Root of Mean Square Error (RMSE): 7.77027214759608
Correlation Coef. (R): 0.999816870896314
R-Square: 0.999633775328896

Parameter	Best Estimate       
---------	-------------       
rsolc    	30.774676380514     
cm       	-4.99365889506396E-5
rq       	-569988.599981876   
rn       	-20379.2697022919   
ln      	8935.93817911551    
rkch     	-10793.0313410194   
rm       	-185.85114389438    
lm       	2.27934259827143    
rl       	167.792508228924    
rctc     	110.815314265582    
t        	1.99088288553436E-13
qc       	37.7983339540601