Why do I get different values when I run gamultiobj for my curve-fitting problem?

Hello Everyone,

I am trying to solve a curve-fitting problem. Previously I was using lsqcurvefit for this task and it worked quite well but since I was worried about local minimum I switched to gamultiobj to avoid this problem. Now I have another problem that I get different values for my B matrix every time that I run gamultiobj (It doesn't give me unique values). What should I do to get unique values for my problem? Are these values still local minimum?Here is my code:

clc
clear all;
close all;

%experimental data (the curve that I am trying to fit to)
[d,s,r] = xlsread('HIT_webplotdigitizer.csv');
t = d(:,1);
y = d(:,2);


%-----------------------------------------------------------------------------------------------------------------------
% Initial guess for B matrix
B0 =[0.8,0.09,0.09,0.09,0.09,0.09,0.09,0.09,0];

m0=size(B0,2);

FitnessFunction = @(B) norm(y-HIT_curve_fitting_4_Gaussian(B, t));
numberOfVariables = m0; % Number of decision variables
lb =[0,0,0,0,0,0,0,0,0];%Lower bound
 ub=[Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf]; % Upper bound
A = []; % No linear inequality constraints
b = []; % No linear inequality constraints
Aeq = []; % No linear equality constraints
beq = []; % No linear equality constraints
  

options = optimoptions(@gamultiobj,'CreationFcn',{@gacreationlinearfeasible},'FunctionTolerance',1e-8,...
         'MaxGenerations',600,'MaxStallGenerations',400,'PopulationSize',700,'MutationFcn',{@mutationadaptfeasible},...
         'CrossoverFraction',0.8,'CrossoverFcn',{@crossoverintermediate});

[B,Fval,exitFlag,Output] = gamultiobj(FitnessFunction,numberOfVariables,A,b,Aeq,beq,lb,ub,options);
% [B,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(odefun1,B0,xx,yy,lb,ub);

figure;
plot(t,y,t,HIT_curve_fitting_4_Gaussian(B, t),'LineWidth',2)
legend('Real HIT','curve-fitting')

function total = HIT_curve_fitting_4_Gaussian(B, x1)
%set temp vector
t0=x1(1);
tf=x1(end);
N=10000;
dt=(tf-t0)/N;
% t=t0:dt:tf;

%----------------------------------------------------------------------------------------------
%function f1
%Second function of the convolution
f1=0.005149*exp(-((x1-62.71)./0.4469).^2);


%construct the firxt part of the convolution
y1=B(1)+B(2)*(x1);


%Convolution of these 2 functions
z1=conv(y1,f1,'full');
con1 = z1(1:length(x1))*dt;
%---------------------------------------------------------------------------------------------
%function f2
%Second function of the convolution
f2=0.01566*exp(-((x1-63.09)./1.209).^2);


%construct the firxt part of the convolution
y2=B(1)+B(3)*(x1);



%Convolution of these 2 functions
z2=conv(y2,f2,'full');
con2 = z2(1:length(x1))*dt;

%----------------------------------------------------------------------------------------
%function f3
%Second function of the convolution
f3=0.02416*exp(-((x1-61.79)./0.8485).^2);



%construct the firxt part of the convolution
y3=B(1)+B(4)*(x1);



%Convolution of these 2 functions
z3=conv(y3,f3,'full');
con3 = z3(1:length(x1))*dt;

%----------------------------------------------------------------------------------------
%function f4
%Second function of the convolution
f4=0.03507*exp(-((x1-63.81)./2.287).^2);

%construct the firxt part of the convolution
y4=B(1)+B(5)*(x1);

%Convolution of these 2 functions
z4=conv(y4,f4,'full');
con4 = z4(1:length(x1))*dt;

%----------------------------------------------------------------------------------------
%function f5
%Second function of the convolution
f5=0.01613*exp(-((x1-70.82)./5.232).^2);



%construct the firxt part of the convolution
y5=B(1)+B(6)*(x1);



%Convolution of these 2 functions
z5=conv(y5,f5,'full');
con5 = z5(1:length(x1))*dt;

%----------------------------------------------------------------------------------------
%function f6
%Second function of the convolution
f6=0.02314*exp(-((x1-66.63)./3.445 ).^2);

%construct the firxt part of the convolution
y6=B(1)+B(7)*(x1);

%Convolution of these 2 functions
z6=conv(y6,f6,'full');
con6 = z6(1:length(x1))*dt;

%----------------------------------------------------------------------------------------
%function f7
%Second function of the convolution
f7=0.001766*exp(-((x1-78.75)./2.143).^2);


%construct the firxt part of the convolution
y7=B(1)+B(8)*(x1);


%Convolution of these 2 functions
z7=conv(y7,f7,'full');
con7 = z7(1:length(x1))*dt;

%----------------------------------------------------------------------------------------
%Create plots
f=f1+f2+f3+f4+f5+f6+f7;
total=con1+con2+con3+con4+con5+con6+con7+B(9);
end

I have attached HIT_webplotdigitizer.csv as well

Answer 1

[Matt J]

I'm not sure why you thought gamultiobj was appropriate for this. gamultiobj is a multi-objective optimizer and therefore will return a range of B values (see also the documentation) across the Pareto front. I suspect you really intended to use ga() if you were looking for a global optimizer.

However, whatever problem made you abaondon lsqcurvefit, I don't think it has to do with local minima. Your HIT_curve_fitting_4_Gaussian looks like a linear function of B, and therefore your overall fitness function should be convex. You could even use lsqlin to find the optimal B, if you could express the linear function HIT_curve_fitting_4_Gaussian in matrix form. That's something you should be able to do quite easily with func2mat on the File Exchange.

https://www.mathworks.com/matlabcentral/fileexchange/44669-func2mat-convert-linear-function-to-matrix