Curve fitting by Genetic Algorithm
Show older comments
Hi everybody, I have a very simple theoric problem but I don't know how to resolve it.
I just wanna do a curve fitting like this:
I want the red line to fit the red dashed one as much as possible in the [0.3 - 0.6] interval.
The red line is a 100x1 vector (x-axis) and 100x1 vector (y-axis).
The red dashed curve is given by this equation (much longer, cropped here):
function [ d2 ] = d2_current( alpha,beta,phi,N,V )
d2 = (4903985730770845.*N.*(((exp(alpha.*(phi - V.*(beta - 1))) + 1).*((alpha.^2.*beta.^2.*exp(alpha.*...
end
Where V is the dependent variable, the previous 100x1 vector (x-axis) corresponding to the [0.3 - 0.6] interval.
The equation depends on 4 independent parameters, each one of them may be included in these intervals:
beta; from to a 1
alpha; from 0 to 15
N; from 0 to 300
phi; from 0 to 3
I have to vary these parameters in order to obtain the best fitting, note that the parameters value must be the same for one single fitting.
The first thing I did was 4 different 'for' loops for each parameter and try to compare each single result with the red curve my means of the euclidean distance. (simplified code):
cont = 0;
for beta = 0:0.01:1
for alpha = 0:0.01:15
for N = 0:1:300
for phi = 0:0.01:3
distance_euclidean(cont,1) = norm(current_red_curve - current_red_dashed_curve);
cont = cont + 1;
end
end
end
end
minim_eucl = min(distance_euclidean);
The minimum euclidean distance would be the best fitting. For the graph example:
beta = 1;
alpha = 15;
N = 300;
phi = 0.45;
Everything is fine until this point, but there's a big problem: the computing time is just huge when the step in the 4 intervals decreases (needed because of the poor fitting).
Looking for better solutions I found the possible way to go: the genetic algorithm. I've been trying it on the Matlab optimization tool, but with no results.
Could somebody help me out?
Thank you very much for your help.
Regards,
8 Comments
Star Strider
on 18 Oct 2019
Post your data, and the correct ‘d2_current’ function.
What are your independent and dependent variables? What parameters do you want to estimate?
Juan Castillo
on 18 Oct 2019
Juan Castillo
on 22 Oct 2019
Update: Little issue here.
Apart from the original fitting, I have to do a second one, but having in mind that the parameters selected must be the optimal for both functions at the same time.
Fitting for my first function:
X = string([0,301067677000000
0,314157576000000
0,327247475000000
0,340337374000000
0,353427273000000
0,366517172000000
0,379607071000000
0,392696970000000
0,405786869000000
0,418876768000000
0,431966667000000
0,445056566000000
0,458146465000000
0,471236364000000
0,484326263000000
0,497416162000000
0,510506061000000
0,523595960000000
0,536685859000000
0,549775758000000
0,562865657000000
0,575955556000000
0,589045455000000
0,602135354000000]);
Y = string([0.0193930890000000
0.0221282684000000
0.0245255339000000
0.0265479770000000
0.0281674962000000
0.0293643273000000
0.0301265743000000
0.0304497403000000
0.0303362581000000
0.0297950211000000
0.0288409140000000
0.0274943436000000
0.0257807698000000
0.0237302359000000
0.0213768998000000
0.0187585647000000
0.0159162099000000
0.0128935212000000
0.00973642228000000
0.00649260510000000
0.00321106079000000
-5.83896230000000e-05
-0.00326556436000000
-0.00636039014000000]);
V = sscanf(sprintf('%s.%s\n',X.'),'%f');
y = str2double(Y);
d2 = @(beta,alpha,N,phi,V) (4903985730770845.*N.*(((exp(alpha.*(phi - V.*(beta - 1))) + 1).*((alpha.^2.*beta.^2.*exp(alpha.*...
(phi - beta.*V)))./(exp(alpha.*(phi - V.*(beta - 1))) + 1) - (alpha.^2.*exp(alpha.*(phi - V.*(beta - 1))).*...
(exp(alpha.*(phi - beta.*V)) + 1).*(beta - 1).^2)./(exp(alpha.*(phi - V.*(beta - 1))) + 1).^2 + ...
(2.*alpha.^2.*exp(2.*alpha.*(phi - V.*(beta - 1))).*(exp(alpha.*(phi - beta.*V)) + 1).*(beta - 1).^2)./...
(exp(alpha.*(phi - V.*(beta - 1))) + 1).^3 - (2.*alpha.^2.*beta.*exp(alpha.*(phi - beta.*V)).*exp(alpha.*...
(phi - V.*(beta - 1))).*(beta - 1))./(exp(alpha.*(phi - V.*(beta - 1))) + 1).^2))./(alpha.*(exp(alpha.*(phi - beta.*V))...
+ 1)) + (exp(alpha.*(phi - V.*(beta - 1))).*(beta - 1).*((alpha.*beta.*exp(alpha.*(phi - beta.*V)))./(exp(alpha.*(phi -...
V.*(beta - 1))) + 1) - (alpha.*exp(alpha.*(phi - V.*(beta - 1))).*(exp(alpha.*(phi - beta.*V)) + 1).*(beta - 1))./...
(exp(alpha.*(phi - V.*(beta - 1))) + 1).^2))./(exp(alpha.*(phi - beta.*V)) + 1) - (beta.*exp(alpha.*(phi - beta.*V)).*...
((alpha.*beta.*exp(alpha.*(phi - beta.*V)))./(exp(alpha.*(phi - V.*(beta - 1))) + 1) - (alpha.*exp(alpha.*(phi - V.*...
(beta - 1))).*(exp(alpha.*(phi - beta.*V)) + 1).*(beta - 1))./(exp(alpha.*(phi - V.*(beta - 1))) + 1).^2).*(exp(alpha.*...
(phi - V.*(beta - 1))) + 1))./(exp(alpha.*(phi - beta.*V)) + 1).^2))./63407003003326144512;
Ve = V(1:18);
ye = y(1:18);
ftns = @(b) norm(ye - d2(b(1),b(2),b(3),b(4),Ve));
PopSz = 500;
Parms = 4;
opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-2, 'MaxGenerations',2E3, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);
t0 = clock;
fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)
[B,fval,exitflag,output] = ga(ftns, Parms, [],[],[],[],[0 0 0 0],[1 15 300 3],[],[],opts)
t1 = clock;
fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)
GA_Time = etime(t1,t0)
%fprintf('\nElapsed Time: %23.15E\t\t%02d:%02d:%02d.%03d\n', GA_Time, hmsdv(GA_Time))
fprintf(1,'\tParameters:\n')
for k1 = 1:length(B)
fprintf(1, '\t\tB(%d) = %8.5f\n', k1, B(k1))
end
xv = linspace(min(Ve), max(Ve), 250);
figure
plot(Ve, ye, 'p')
hold on
plot(xv, d2(B(1),B(2),B(3),B(4),xv), '-r')
grid
Perfect until here, but I wanna do another fitting at the same time, so I add the new function and its Vc and Ic values as X and Y axis where I want the fitting:
X = string([0,301067677000000
0,314157576000000
0,327247475000000
0,340337374000000
0,353427273000000
0,366517172000000
0,379607071000000
0,392696970000000
0,405786869000000
0,418876768000000
0,431966667000000
0,445056566000000
0,458146465000000
0,471236364000000
0,484326263000000
0,497416162000000
0,510506061000000
0,523595960000000
0,536685859000000
0,549775758000000
0,562865657000000
0,575955556000000
0,589045455000000
0,602135354000000]);
Y = string([0.0193930890000000
0.0221282684000000
0.0245255339000000
0.0265479770000000
0.0281674962000000
0.0293643273000000
0.0301265743000000
0.0304497403000000
0.0303362581000000
0.0297950211000000
0.0288409140000000
0.0274943436000000
0.0257807698000000
0.0237302359000000
0.0213768998000000
0.0187585647000000
0.0159162099000000
0.0128935212000000
0.00973642228000000
0.00649260510000000
0.00321106079000000
-5.83896230000000e-05
-0.00326556436000000
-0.00636039014000000]);
V = sscanf(sprintf('%s.%s\n',X.'),'%f');
y = str2double(Y);
d2 = @(beta,alpha,N,phi,V) (4903985730770845.*N.*(((exp(alpha.*(phi - V.*(beta - 1))) + 1).*((alpha.^2.*beta.^2.*exp(alpha.*...
(phi - beta.*V)))./(exp(alpha.*(phi - V.*(beta - 1))) + 1) - (alpha.^2.*exp(alpha.*(phi - V.*(beta - 1))).*...
(exp(alpha.*(phi - beta.*V)) + 1).*(beta - 1).^2)./(exp(alpha.*(phi - V.*(beta - 1))) + 1).^2 + ...
(2.*alpha.^2.*exp(2.*alpha.*(phi - V.*(beta - 1))).*(exp(alpha.*(phi - beta.*V)) + 1).*(beta - 1).^2)./...
(exp(alpha.*(phi - V.*(beta - 1))) + 1).^3 - (2.*alpha.^2.*beta.*exp(alpha.*(phi - beta.*V)).*exp(alpha.*...
(phi - V.*(beta - 1))).*(beta - 1))./(exp(alpha.*(phi - V.*(beta - 1))) + 1).^2))./(alpha.*(exp(alpha.*(phi - beta.*V))...
+ 1)) + (exp(alpha.*(phi - V.*(beta - 1))).*(beta - 1).*((alpha.*beta.*exp(alpha.*(phi - beta.*V)))./(exp(alpha.*(phi -...
V.*(beta - 1))) + 1) - (alpha.*exp(alpha.*(phi - V.*(beta - 1))).*(exp(alpha.*(phi - beta.*V)) + 1).*(beta - 1))./...
(exp(alpha.*(phi - V.*(beta - 1))) + 1).^2))./(exp(alpha.*(phi - beta.*V)) + 1) - (beta.*exp(alpha.*(phi - beta.*V)).*...
((alpha.*beta.*exp(alpha.*(phi - beta.*V)))./(exp(alpha.*(phi - V.*(beta - 1))) + 1) - (alpha.*exp(alpha.*(phi - V.*...
(beta - 1))).*(exp(alpha.*(phi - beta.*V)) + 1).*(beta - 1))./(exp(alpha.*(phi - V.*(beta - 1))) + 1).^2).*(exp(alpha.*...
(phi - V.*(beta - 1))) + 1))./(exp(alpha.*(phi - beta.*V)) + 1).^2))./63407003003326144512;
Ve = V(1:18);
ye = y(1:18);
Vc = [0;0.0130898990000000;0.0261797980000000;0.0392696970000000;0.0523595960000000;0.0654494949000000;0.0785393939000000;0.0916292929000000;0.104719192000000;0.117809091000000;0.130898990000000;0.143988889000000;0.157078788000000;0.170168687000000;0.183258586000000;0.196348485000000;0.209438384000000;0.222528283000000;0.235618182000000;0.248708081000000;0.261797980000000;0.274887879000000;0.287977778000000;0.301067677000000;0.314157576000000;0.327247475000000;0.340337374000000;0.353427273000000;0.366517172000000;0.379607071000000;0.392696970000000];
Ic = [4.53936181000000e-05;0.000399955721000000;0.000770116251000000;0.00115194119000000;0.00154212253000000;0.00193791839000000;0.00233709646000000;0.00273788056000000;0.00313890039000000;0.00353914431000000;0.00393791499000000;0.00433478809000000;0.00472957355000000;0.00512227977000000;0.00551308032000000;0.00590228321000000;0.00629030269000000;0.00667763337000000;0.00706482673000000;0.00745246979000000;0.00784116597000000;0.00823151805000000;0.00862411304000000;0.00901950904000000;0.00941822390000000;0.00982072562000000;0.0102274244000000;0.0106386664000000;0.0110547288000000;0.0114758166000000;0.0119020602000000];
curr = @(beta,alpha,N,phi,V) 2*(1.602176462e-19)^2*N/(6.62606876e-34)*(V+1/alpha*log((1+exp(alpha*(phi-beta*V)))/(1+exp(alpha*(phi+(1-beta)*V)))));
% Minimizing the distance of both functions at the same time
ftns = @(b) norm(ye - d2(b(1),b(2),b(3),b(4),Ve) + norm(Ic - curr(b(1),b(2),b(3),b(4),Vc)));
PopSz = 500;
Parms = 4;
opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-2, 'MaxGenerations',2E3, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);
t0 = clock;
fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)
[B,fval,exitflag,output] = ga(ftns, Parms, [],[],[],[],[0 0 0 0],[1 15 300 3],[],[],opts)
t1 = clock;
fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)
GA_Time = etime(t1,t0)
%fprintf('\nElapsed Time: %23.15E\t\t%02d:%02d:%02d.%03d\n', GA_Time, hmsdv(GA_Time))
fprintf(1,'\tParameters:\n')
for k1 = 1:length(B)
fprintf(1, '\t\tB(%d) = %8.5f\n', k1, B(k1))
end
xv = linspace(min(Ve), max(Ve), 250);
figure
plot(Ve, ye, 'p')
hold on
plot(xv, d2(B(1),B(2),B(3),B(4),xv), '-r')
grid
This is what I get:
The new fitting is 0. Ok, no problem, let's fit only the new function just to try. So the line now looks like this:
ftns = @(b) norm(Ic - curr(b(1),b(2),b(3),b(4),Vc));
And adding the new plot:
plot(Vc, Ic, '-k')
This is the result:
Still 0. Here the red line should be along with the black one (from 0.3 to 0.4).
The main goal here is look for the parameters which best fit both the black and the blue one. What am I doing wrong?
Many thanks in advance.
Star Strider
on 22 Oct 2019
‘The main goal here is look for the parameters which best fit both the black and the blue one.’
I doubt that is possible. The (Ve,ye) data plot a curve (blue pentagrams in your plot), and the (Vc,Ic) data plot a straight line (the black line in your plot).
‘What am I doing wrong?’
You are assuming that the parameters for your function that describe the blue pentagram data curve will also fit the black data straight line. They have absolutely nothing in common. If you have two data sets that were similar, with similar ranges and other characteristics in both the independent and dependent variables, you could probably get a parameter set that fit both of them (although neither as well as specific parameter estimates for each one). However here that is not possible.
Juan Castillo
on 23 Oct 2019
Star Strider
on 23 Oct 2019
The expression itself appears to me to be correct, although I have never done such an optimization. The problem of course is that you are attempthing to identify the same parameters from two different data sets of different lengths, with different characteristics, and over different ranges of the indpendent variable. I cannot imagine that this will ever be successful.
Juan Castillo
on 23 Oct 2019
Star Strider
on 23 Oct 2019
As always, my pleasure.
Accepted Answer
Star Strider
on 18 Oct 2019
I can fit part of your data, however not all of it, because there is some noise in elements 19 through 24 that your function is likely to not able to account for. (I had to eliminate the commas from your data and translate them into decimal points. Thjat worked for the data you posted, however I do not know if it will work for your entire vector.)
Try this:
X = string([0,301067677000000
0,314157576000000
0,327247475000000
0,340337374000000
0,353427273000000
0,366517172000000
0,379607071000000
0,392696970000000
0,405786869000000
0,418876768000000
0,431966667000000
0,445056566000000
0,458146465000000
0,471236364000000
0,484326263000000
0,497416162000000
0,510506061000000
0,523595960000000
0,536685859000000
0,549775758000000
0,562865657000000
0,575955556000000
0,589045455000000
0,602135354000000]);
Y = string([0,0193930890000000
0,0221282684000000
0,0245255339000000
0,0265479770000000
0,0281674962000000
0,0293643273000000
0,0301265743000000
0,0304497403000000
0,0303362581000000
0,0297950211000000
0,0288409140000000
0,0274943436000000
0,0257807698000000
0,0237302359000000
0,0213768998000000
0,0187585647000000
0,0159162099000000
0,0128935212000000
0,00973642228000000
0,00649260510000000
0,00321106079000000
-5,83896230000000e-05
-0,00326556436000000
-0,00636039014000000]);
V = sscanf(sprintf('%s.%s\n',X.'),'%f');
y = sscanf(sprintf('%s.%s\n',Y.'),'%f');
d2 = @(beta,alpha,N,phi,V) (4903985730770845.*N.*(((exp(alpha.*(phi - V.*(beta - 1))) + 1).*((alpha.^2.*beta.^2.*exp(alpha.*...
(phi - beta.*V)))./(exp(alpha.*(phi - V.*(beta - 1))) + 1) - (alpha.^2.*exp(alpha.*(phi - V.*(beta - 1))).*...
(exp(alpha.*(phi - beta.*V)) + 1).*(beta - 1).^2)./(exp(alpha.*(phi - V.*(beta - 1))) + 1).^2 + ...
(2.*alpha.^2.*exp(2.*alpha.*(phi - V.*(beta - 1))).*(exp(alpha.*(phi - beta.*V)) + 1).*(beta - 1).^2)./...
(exp(alpha.*(phi - V.*(beta - 1))) + 1).^3 - (2.*alpha.^2.*beta.*exp(alpha.*(phi - beta.*V)).*exp(alpha.*...
(phi - V.*(beta - 1))).*(beta - 1))./(exp(alpha.*(phi - V.*(beta - 1))) + 1).^2))./(alpha.*(exp(alpha.*(phi - beta.*V))...
+ 1)) + (exp(alpha.*(phi - V.*(beta - 1))).*(beta - 1).*((alpha.*beta.*exp(alpha.*(phi - beta.*V)))./(exp(alpha.*(phi -...
V.*(beta - 1))) + 1) - (alpha.*exp(alpha.*(phi - V.*(beta - 1))).*(exp(alpha.*(phi - beta.*V)) + 1).*(beta - 1))./...
(exp(alpha.*(phi - V.*(beta - 1))) + 1).^2))./(exp(alpha.*(phi - beta.*V)) + 1) - (beta.*exp(alpha.*(phi - beta.*V)).*...
((alpha.*beta.*exp(alpha.*(phi - beta.*V)))./(exp(alpha.*(phi - V.*(beta - 1))) + 1) - (alpha.*exp(alpha.*(phi - V.*...
(beta - 1))).*(exp(alpha.*(phi - beta.*V)) + 1).*(beta - 1))./(exp(alpha.*(phi - V.*(beta - 1))) + 1).^2).*(exp(alpha.*...
(phi - V.*(beta - 1))) + 1))./(exp(alpha.*(phi - beta.*V)) + 1).^2))./63407003003326144512;
Ve = V(1:18);
ye = y(1:18);
ftns = @(b) norm(ye - d2(b(1),b(2),b(3),b(4),Ve));
PopSz = 500;
Parms = 4;
opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-2, 'MaxGenerations',2E3, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);
t0 = clock;
fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)
[B,fval,exitflag,output] = ga(ftns, Parms, [],[],[],[],[],[],[],[],opts)
t1 = clock;
fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)
GA_Time = etime(t1,t0)
fprintf('\nElapsed Time: %23.15E\t\t%02d:%02d:%02d.%03d\n', GA_Time, hmsdv(GA_Time))
fprintf(1,'\tParameters:\n')
for k1 = 1:length(B)
fprintf(1, '\t\tB(%d) = %8.5f\n', k1, B(k1))
end
xv = linspace(min(Ve), max(Ve), 250);
figure
plot(Ve, ye, 'p')
hold on
plot(xv, d2(B(1),B(2),B(3),B(4),xv), '-r')
grid
producing these parameter estimates:
Parameters:
B(1) = 15.96980
B(2) = 0.63200
B(3) = 542.18665
B(4) = 7.56829
with the elements of ‘B’ being, in order: ‘beta’, ‘alpha’, ‘N’, and ‘phi’. They retain their full precision internally.
and this plot:
Plot images are not showing up with the ‘Picture’ icon for some reason, so I am also attaching the plot image as a separate file to be sure it is visible.
7 Comments
Juan Castillo
on 18 Oct 2019
Star Strider
on 18 Oct 2019
My pleasure.
You can implement the constraints in ga as:
[B,fval,exitflag,output] = ga(ftns, Parms, [],[],[],[],[0 0 0 0],[1 15 300 3],[],[],opts)
I do not understand:
beta; from [0 to a 1]
since I do not see ‘a’ anywhere. (I just use 0 and 1 here.)
The constraints prevent your function from fitting your data. The best fit parameters (several ga runs) are:
Parameters:
B(1) = 1.00000
B(2) = 15.00000
B(3) = 300.00000
B(4) = 0.40335
The ‘hmsdv’ function is one I wrote about 25 years ago to parse the etime output, long before tic, toc, and other time functions existed. Just eliminate this line that calls it:
fprintf('\nElapsed Time: %23.15E\t\t%02d:%02d:%02d.%03d\n', GA_Time, hmsdv(GA_Time))
or use tic and toc in place of the clock calls.
It is not unusual to get widely-varying parameter values that nevertheless produce a successful approximation of the objective function to the data with ga (and other optimisation algorithms), since different combinations of values will still produce values of the objective function (especially with complicated objective functions) that create approximately the same values within the function. If the fit is appropriate (looking at the plot or using statistical goodness-of-fit measures), then the parameters are likely correct.
Since the parameter constraints appear to be necessary in order to be appropriate to the physics of the ‘d2’ function, look carefully at your data to be certain the units, scale factors, and other necessary values match those that ‘d2’ uses and returns. That is the only reason I can think of for the fit to be as poor as it is.
Juan Castillo
on 18 Oct 2019
Star Strider
on 18 Oct 2019
As always, my pleasure!
I very much appreciate your compliments!
Star Strider
on 19 Oct 2019
I did not realise before that my attempt to change the (,) to (.) in your data did not work correctly, (the values were off by an order of magnitude), so I changed them manually.
I then re-ran my ga code:
[B,fval,exitflag,output] = ga(ftns, Parms, [],[],[],[],[0 0 0 0],[1 15 300 3],[],[],opts)
and got these parameters:
Parameters:
B(1) = 0.93329
B(2) = 15.00000
B(3) = 120.59092
B(4) = 0.36837
and this plot:
Juan Castillo
on 19 Oct 2019
Star Strider
on 19 Oct 2019
As always, my pleasure!
More Answers (1)
Alex Sha
on 19 Oct 2019
1 vote
Rerfer the results below:
Root of Mean Square Error (RMSE): 0.00353350868235726
Sum of Squared Residual: 0.00029965640659906
Correlation Coef. (R): 0.975244832242099
R-Square: 0.95110248281492
Adjusted R-Square: 0.946445576416342
Determination Coef. (DC): 0.902466880697908
Chi-Square: -0.0236084937698218
F-Statistic: 41.4027542316901
Parameter Best Estimate
---------- -------------
beta1 1
alpha 15
n 105.578188880725
phi 0.388542678989983
1 Comment
Juan Castillo
on 19 Oct 2019
Categories
Find more on Solver Outputs and Iterative Display in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
