Non linear fitting and parametric optimization with genetic algorithm

I have a constitutive model with 4 parameters to fit to experimental stress vs strain data. Before implementing the GA, I tried lsqnonlin (with and without multistart), but one of the main problems is that the algorithm often tends to assign bound values to the optimised parameters. To the first parameter, which is larger than the others, it assigns the upper bound, to the second and third it assigns the lower bounds. Plotting the best individuals over time, it appears that the ga spans little between the solutions, and the optimisation seems to be biased towards the first parameter.
I am using the uniaxial model proposed by Dong, and have rechecked the correctness of the model form several times. Below is the main containing the ga and my error function. Any help is welcome, thank you

%CHANGE LINE 3 TO SELECT THE PATCH
%CHANGE LINES 9 AND 10 TO SELECT THE SPECIMEN
clc; clear all;close all;
load ("data_11122.mat"); % change lot number if necessary
clearvars -except data

%Model requires lambda= x/x0 instead of epsilon: calculate lambda
% from exp data and store in data struct

stress=data(6).stress; %change index to change dogbone if necessary
stress=stress/max(stress);
strain=data(6).strain; %change index to change dogbone if necessary
p_initial=[50 25 5 .5]*10e-3;
a0=[1 0 0];
%vincoli su zeta imposti da dong ( zeta tra 0 e 1 )
lb=[0 0 0 0]; %% AGGIUSTA UNITA DIN MISURA
ub=[650 98 35 1];

%ottimizzazione non lineare con metodo least square ( like Dong )
options = optimoptions(‘ga’, ‘Display’, ‘iter’, ‘PopulationSize’, 100, ‘MaxGenerations’, 200, ‘UseParallel’, true,’NonlinearConstraintAlgorithm’,’penalty’,’PlotFcn’,’gaplotbestindiv’,’InitialPopulationMatrix’,p_initial);
nvars=4;
[p_opt, fval, exitflag, output] = ga(@(p)(error_function(p, strain, stress)), nvars, [], [], [], [], lb, ub, [], options);

%check plot e ricostruzione del modello con i parametri ottimizzati
c = p_opt(1); % c
k1 = p_opt(2); % k1
k2 = p_opt(3); % k2
zeta = p_opt(4); % zeta
for i=1:length(strain)

sigma_calculated(i)=c*(strain(i)^2) + 2*strain(i)^2*k1*((strain(i)^2-1)*(i-zeta)^2*exp(k2*(strain(i)^2-1)^2*(1-zeta)^2));

end

% for i=1:length(strain)
% f = [strain(i) 0 0; %deformation gradient for each value of strain
% 0 1/sqrt(strain(i)) 0;
% 0 0 1/sqrt(strain(i))];
% C = f .* f’; % Calcola il tensore di Cauchy-Green
% prodot = a0 .* a0′;
% I4 = sum(sum(C .* prodot)); % calcolo invariante di C
% I = eye(3); % Matrice identità
%
% sigma= f*c*I*f’ + 2*f*(k1 * (1 – zeta)^2 * (I4 – 1) * (exp(k2 * ((1 – zeta) * (I4 – 1))^2)) * prodot)*f’;
%
% sigma_calculated(i)= sigma(1,1); %componente xx
% end

figure()
plot (strain,stress);
hold on
plot(strain, sigma_calculated,’r’);
legend (‘experimental’, ‘model’,’Location’, ‘northwest’);

function error = error_function(p,strain,stress)

sigma_calculated = zeros(length(strain), 1);

for i = 2:length(strain) %avoiding first value=0;

% f = [strain(i) 0 0; %deformation gradient
% 0 1/sqrt(strain(i)) 0;
% 0 0 1/sqrt(strain(i))];
% C = f .* f’; % tensore di Cauchy-Green
% prodot = a0 .* a0′;
% I4 = a0′ .* C .* a0; % calcolo IV invariante di C
% I = eye(3);
%
% sigma= f*p(1)*I*f’ + 2*f*(p(2) * (1 – p(4))^2 * (I4 – 1) * (exp(p(3) * ((1 – p(4)) * (I4 – 1))^2)) * prodot)*f’;
%
%
% sigma_calculated(i)= sigma(1,1); %componente xx
sigma_calculated(i)=p(1)*(strain(i)^2) + 2*strain(i)^2*p(2)*((strain(i)^2-1)*(i-p(4))^2*exp(p(3)*(strain(i)^2-1)^2*(1-p(4))^2));

end

error= sum((sigma_calculated-stress).^2) ; %/sum(stress).^2; % Somma dei quadrati delle differenze
%aggiungere normalizzazione rispetto allo stress( guardare script marta calcolo error )
endI have a constitutive model with 4 parameters to fit to experimental stress vs strain data. Before implementing the GA, I tried lsqnonlin (with and without multistart), but one of the main problems is that the algorithm often tends to assign bound values to the optimised parameters. To the first parameter, which is larger than the others, it assigns the upper bound, to the second and third it assigns the lower bounds. Plotting the best individuals over time, it appears that the ga spans little between the solutions, and the optimisation seems to be biased towards the first parameter.
I am using the uniaxial model proposed by Dong, and have rechecked the correctness of the model form several times. Below is the main containing the ga and my error function. Any help is welcome, thank you

%CHANGE LINE 3 TO SELECT THE PATCH
%CHANGE LINES 9 AND 10 TO SELECT THE SPECIMEN
clc; clear all;close all;
load ("data_11122.mat"); % change lot number if necessary
clearvars -except data

%Model requires lambda= x/x0 instead of epsilon: calculate lambda
% from exp data and store in data struct

%check plot e ricostruzione del modello con i parametri ottimizzati
c = p_opt(1); % c
k1 = p_opt(2); % k1
k2 = p_opt(3); % k2
zeta = p_opt(4); % zeta
for i=1:length(strain)

sigma_calculated(i)=c*(strain(i)^2) + 2*strain(i)^2*k1*((strain(i)^2-1)*(i-zeta)^2*exp(k2*(strain(i)^2-1)^2*(1-zeta)^2));

end

figure()
plot (strain,stress);
hold on
plot(strain, sigma_calculated,’r’);
legend (‘experimental’, ‘model’,’Location’, ‘northwest’);

function error = error_function(p,strain,stress)

sigma_calculated = zeros(length(strain), 1);

for i = 2:length(strain) %avoiding first value=0;

end

error= sum((sigma_calculated-stress).^2) ; %/sum(stress).^2; % Somma dei quadrati delle differenze
%aggiungere normalizzazione rispetto allo stress( guardare script marta calcolo error )
end I have a constitutive model with 4 parameters to fit to experimental stress vs strain data. Before implementing the GA, I tried lsqnonlin (with and without multistart), but one of the main problems is that the algorithm often tends to assign bound values to the optimised parameters. To the first parameter, which is larger than the others, it assigns the upper bound, to the second and third it assigns the lower bounds. Plotting the best individuals over time, it appears that the ga spans little between the solutions, and the optimisation seems to be biased towards the first parameter.
I am using the uniaxial model proposed by Dong, and have rechecked the correctness of the model form several times. Below is the main containing the ga and my error function. Any help is welcome, thank you

%CHANGE LINE 3 TO SELECT THE PATCH
%CHANGE LINES 9 AND 10 TO SELECT THE SPECIMEN
clc; clear all;close all;
load ("data_11122.mat"); % change lot number if necessary
clearvars -except data

%Model requires lambda= x/x0 instead of epsilon: calculate lambda
% from exp data and store in data struct

%check plot e ricostruzione del modello con i parametri ottimizzati
c = p_opt(1); % c
k1 = p_opt(2); % k1
k2 = p_opt(3); % k2
zeta = p_opt(4); % zeta
for i=1:length(strain)

sigma_calculated(i)=c*(strain(i)^2) + 2*strain(i)^2*k1*((strain(i)^2-1)*(i-zeta)^2*exp(k2*(strain(i)^2-1)^2*(1-zeta)^2));

end

figure()
plot (strain,stress);
hold on
plot(strain, sigma_calculated,’r’);
legend (‘experimental’, ‘model’,’Location’, ‘northwest’);

function error = error_function(p,strain,stress)

sigma_calculated = zeros(length(strain), 1);

for i = 2:length(strain) %avoiding first value=0;

end

Cart

Cart