Nonlinear curve fitting problem with parameters to be optimised

in a non-linear fitting context, the output parameters of an optimisation algorithm assume upper bounds and lower bounds, making the fitting unsatisfactory.

In particular, the first parameter takes on the value of upper bounds, while the remaining parameters take on the values of lower bounds.

Below is the main, the error function and a graph showing the unsatisfactory fitting.

Main:

%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(1).stress; %change index to change dogbone if necessary
strain=data(1).strain; %change index to change dogbone if necessary
p_initial=[0.5 0.5 0.5 0.5];
a0=[1 0 0];
%vincoli su zeta imposti da dong ( zeta tra 0 e 1 )
lb=[0.1 20 0.5 0.1];
ub=[150 150 100 1];

%ottimizzazione non lineare con metodo least square ( like Dong )
options = optimoptions('lsqnonlin', 'Display', 'iter', 'Algorithm','levenberg-marquardt');
[p_opt, resnorm,output] = lsqnonlin(@(p) error_function(p,strain,stress,a0), p_initial,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)
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 = 2 * f * (c / 2 * I + 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(1)
plot (strain,stress);
hold on
plot(strain, sigma_calculated,'r');
legend ('experimental', 'model','Location', 'northwest');

%error function

function error = error_function(p,strain,stress,a0)

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 = 2 * f * (p(1) / 2 * I + 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

end

error= sum((sigma_calculated - stress).^2); % Somma dei quadrati delle differenze

end