Nonlinear least-squares data fit

I am trying to make a data fit for the data attached to this post,Nu=f(Re,Theta,Beta).I use lsqnonlin(fun,x0) function for this purpose.I have created a script file for this fitting,but everytime I try to run the script,the program always shows error messages.So,what is the problem with this script.

clc
clear all
% Create an anonymous function that describes the expected relationship
% between X and Y
f=@(c,x) c(1).*(x(:,1).^c(2)).*(x(:,2).^c(3)).*(x(:,3).^c(4))./x(:,4)-1;
% data set 
% Specify x variables from data file,Re,Theta and Beta columns.
x=xlsread('all data for fitting');
% Specify y variable from data file ,(Nu)column.
y=x(:,4);
% Specify a vector of starting conditions for the solvers
c0=[1;1;1;1];
% Perform a nonlinear regression
c=lsqnonlin(f,c0);

Answer 1

[Star Strider]

The objective function needs to be coded as:

ffcn = @(c) f(c,x) - y; 

with the complete lsqnonlin call being:

f=@(c,x) c(1).*(x(:,1).^c(2)).*(x(:,2).^c(3)).*(x(:,3).^c(4))./x(:,4)-1;
ffcn = @(c) f(c,x) - y; 
c0=[1;1;1;1];
C = lsqnonlin(ffcn, c0);

producing:

C =
   1.0308e-01
   1.3246e+00
   1.9801e-06
  -4.6017e-01

.

[MOHAMED ABDULAZIM]

I have already made this step in defining the function,I aimed to form the function in the relative error form ((Fit -actual) /actual).

[Star Strider]

That was not part of the Question you asked.

I do not understand what you want to do.

I used the function you provided. How does that differ from what you want?

This version of the objecctive function runs:

ffcn = @(c) (f(c,x) - y)./y; 

and produces approximately the same parameter values.

[MOHAMED ABDULAZIM]

you are correct about the objective function and about what you have said,but I mean that the function that you wrote (ffcn = @(c) (f(c,x) - y)./y;),I have already created it in the code(f=@(c,x) c(1).*(x(:,1).^c(2)).*(x(:,2).^c(3)).*(x(:,3).^c(4))./x(:,4)-1) ,the form is different as I simplified the fraction.What is the error that make the code stop from running?

[MOHAMED ABDULAZIM]

I tried the code with the modifications that you have suggested but there is a message after running the script.

Local minimum possible.

lsqnonlin stopped because the size of the current step is less than
the default value of the step size tolerance.

<stopping criteria details>

clc
clear all
% Create an anonymous function that describes the expected relationship
% between X and Y
f=@(c,x) c(1).*(x(:,1).^c(2)).*(x(:,2).^c(3)).*(x(:,3).^c(4));
% data set 
% Specify x variables from data file,Re,Theta and Beta columns.
x=xlsread('all data for fitting');
% Specify y variable from data file ,(Nu)column.
y=x(:,4);
% objective function 
ff=@(c)(f(c,x)-y)./y;
% Specify a vector of starting conditions for the solvers
c0=[1;1;1;1];
% Perform a nonlinear regression
c=lsqnonlin(ff,c0);

Here is the code with the modifications

[Star Strider]

That is perfectly normal. The lsqnonliln message is not an error, simply letting you know that while it cocnverged, it may not have converged on the global minimum.

I ran this with the genetic algorithm function seveeral times, and the best parameter set it found was:

4.8773e+000   479.3000e-003    16.3000e-003  -555.5538e-003

with a fitness of 5.5151. That is likely as close as it is possible to get to the global minimum.

The lsqnonlin function improved on those slightly to return:

4.3202e+000   496.4773e-003    21.7014e-006  -543.8021e-003

Experiment with those values for the parameters.

.

[MOHAMED ABDULAZIM]

How to use the genetic algorithm function in this fitting?

[Star Strider]

That is a bit difficult to describe.

I used this:

D = xlsread('all data for fitting.xlsx');
x = D;
y = x(:,4);
f=@(c,x) c(1).*(x(:,1).^c(2)).*(x(:,2).^c(3)).*(x(:,3).^c(4));
ffcn = @(c) (f(c,x) - y)./y; 
ftns = @(c) norm(ffcn(c));

PopSz = 500;
Parms = 4;
opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-4, 'MaxGenerations',2E3, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);
t0 = clock;
fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)
[theta,fval,exitflag,output] = ga(ftns, Parms, [],[],[],[],-Inf(Parms,1),Inf(Parms,1),[],[],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 s', GA_Time)

fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%d) = %12.5E\n', k1, theta(k1))
end

The ga function requires the Global Optimization Toolbox. This code should run in the R2013a version, however I cannot test that, so I cannot be certain.

[MOHAMED ABDULAZIM]

Did you use the genetic algorithm and lsqnonlin function separately for comparison purpose?

[Star Strider]

I used the genetic algorithm ga function first, then used its returned parameter vector as the initial parameter estimates for lsqnonlin.

[MOHAMED ABDULAZIM]

I have tried to take the initial values from genetic algorithm and use them in the least square nonlinear fit but every time I try to run the script,an error message appears.

clc
clear all
% Create an anonymous function that describes the expected relationship
% between X and Y
f=@(c,x) c(1).*(x(:,1).^c(2)).*exp((x(:,3)*c(4))+(x(:,2)*c(3)));
% data set 
% Specify x variables from data file,Re,Theta and Beta columns.
x=xlsread('all data for fitting');
% Specify y variable from data file ,(Nu)column.
y=x(:,4);
% objective function 
ff=@(c)(f(c,x)-y)./y;
% maximum of objective function
ffmax=@(c)norm(ff(c));
% Identifying population size and number of parameters for genetic
% algorithm
PopSz = 500;
Parms = 4;
% Modifying the options of genetic algotithm optimization function
options = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-4, 'MaxGenerations',2E3, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);
% Genetic algorithm application
[theta,fval,exitflag,output] = ga(ffmax, Parms, [],[],[],[],-Inf(Parms,1),Inf(Parms,1),[],[],options);
% Specify a vector of starting conditions for the solvers(Least squre
% nnonlinear fit)
c0=[theta(1) ;theta(2);theta(3); theta(4)];
% Perform a nonlinear regression
c=lsqnonlin(ff,c0);
y1=f(c,x);
diff_Nu=y-y1;
Abs_diff_Nu=abs(diff_Nu);
perc_error=(Abs_diff_Nu./y);
x_line45=0:max(y);y_line45=0:max(y);
[max_error,I_max]=max(perc_error);
[min_error,I_min]=min(perc_error);
x_line_high=0:max(y); y_line_high=(max_error+1)*(0:max(y));
x_line_low=0:max(y); y_line_low=(1-max_error)*(0:max(y));
figure
plot(x_line45,y_line45,'k')
hold on
scatter(y,y1)
hold on
plot(x_line_high,y_line_high,'r')
hold on
plot(x_line_low,y_line_low,'r')
N=[y,y1];
header=['Nuact','Nufit'];
disp(header)
disp(N)

Here is the error:

Error using optimoptions (line 114)
Invalid solver specified. Provide a solver name or handle
(such as 'fmincon' or @fminunc).
Type DOC OPTIMOPTIONS for a list of solvers.

Error in bestfitever (line 20)
options = optimoptions('ga', 'PopulationSize',PopSz,
'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-4,
'MaxGenerations',2E3, 'PlotFcn',@gaplotbestf,
'PlotInterval',1);

[Star Strider]

I have no idea what is throwing that error. The code appears to be correct, however my code runs correctly on R2020a, so be certain the code is correct for R2013a. There have been many changes in the last 7 years, and I no longer have access to the documentation for R2013a.

[MOHAMED ABDULAZIM]

Did you try to run this script on your version?

[Star Strider]

I used this:

D = xlsread('all data for fitting.xlsx');
x = D;
y = x(:,4);
f=@(c,x) c(1).*(x(:,1).^c(2)).*(x(:,2).^c(3)).*(x(:,3).^c(4));
ffcn = @(c) (f(c,x) - y)./y; 
ftns = @(c) norm(ffcn(c));

PopSz = 500;
Parms = 4;
opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-4, 'MaxGenerations',2E3, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);
t0 = clock;
fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)
[theta,fval,exitflag,output] = ga(ftns, Parms, [],[],[],[],-Inf(Parms,1),Inf(Parms,1),[],[],opts)
t1 = clock;
fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)
GA_Time = etime(t1,t0)
QQQ1 = datetime([zeros(1,5) GA_Time], 'Format','HH:mm:ss.SSS')
fprintf('\nElapsed Time: %23.15E s ', GA_Time)

fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%d) = %12.5E\n', k1, theta(k1))
end

and when I ran that just now, got these parameter estimates:

theta =
     2.8517e+000   431.7000e-003    99.6000e-003  -324.6437e-003

with a fitness value of:

fval =
     5.5386e+000

.