Problem with lsqnonlin and error function implementing a bovine pericardium constitutive model

I implemented a constitutive model for the bovine pericardium to be fit with experimental data to derive four model parameters. The code runs, but I have doubts about the correct implementation of the error function. Showing the iterations on the screen, I notice that at certain points the objective function returns nan or inf. I have checked the experimental data and in some cases there are values close to zero, but I don't think this justifies the fact that the function converges, thus giving me the 4 parameters I am looking for, but without obtaining a decent fitting with the experimental data. Moreover, optimized parameters number 2 and 3 are assigned to lower bounds, while the first optimized parameter is very high. Below is the main and the error function.

Main:

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
% stress=stress/max(abs(stress));
% strain=strain/max(abs(strain));
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 0.1 0.1 0.1];
ub=[inf inf inf inf];

%ottimizzazione non lineare con metodo least square ( like Dong )
options = optimoptions('lsqnonlin', 'Display', 'iter', 'Algorithm','levenberg-marquardt');
[p_opt, resnorm] = 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
prod = a0 .* a0';  
I4 = sum(sum(C .* prod)); % 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))) * prod .* 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,k)
   
    sigma_calculated = zeros(length(strain), 1); 
  
    for i = 1:length(strain)
    
    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
    prod = 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)) * prod) .* f';
    sigma_calculated(i)= sigma(1,1);                                           %componente xx

     
    end
%il calcolo dell'errore va dentro o fuori dal for?
   error= sum((sigma_calculated(i) - stress(i))^2);                           % Somma dei quadrati delle differenze 
%   disp(['Iterazione: ', num2str(k), ' - Errore totale attuale: ', num2str(error)]);

end

[Torsten]

Use

error= sum((sigma_calculated - stress).^2); 

instead of

error= sum((sigma_calculated(i) - stress(i))^2); 

And rename your variable "prod". The name should remain reserved for the MATLAB function "prod".

[Nicolò]

Thank you very much for your time. I have followed your advice, but now the code does not run and I get the following error message:

"Error using lsqncommon (line 15)

Objective function is returning undefined values at initial point. lsqnonlin cannot continue.

Error in lsqnonlin (line 260)

lsqncommon(funfcn,xCurrent,lb,ub,options,defaultopt,optimgetFlag,caller,...

Error in main (line 22)

[p_opt, resnorm,output] = lsqnonlin(@(p) error_function(p,strain,stress,a0), p_initial,lb,ub,options); "

Do you have any idea?

[Nicolò]

Edit: i solved the previous messagge error, code runs but i haven't a good fit. The problem on high parameter 1 and low values of parameters 2 and 3 persists. In particular, p(1) takes the ub values, p(2) and p(3) take lb values

[Star Strider]

If you have the Global Optimization Toolbox, experiment by using ga or one of the multiple start functions (MultiStart, GlobalSearch). It may be that your objective function is very sensitive to the initial start point, and your response surface could have multiple local minima.

[Nicolò]

@Star Strider I implemented the multistart function but the result is the same. I notice that the algorithm assigns the upper bound to the first parameter, while it assigns the lower bounds to parameters 2 and 3.

Answer 1

[William Rose]

@Nicolò, can you provide file data_11122.mat?

[Nicolò]

Here you are. note that the 0 values in the dataset are already modified ( it didn't worked although this adjustement)

Answer 2

[Torsten]

strain(1) = 0 - thus you must remove the first data point since you divide by strain(i).

Maybe your model function does not allow a curve with turning point - you should check this for some choices of the parameters. If it does, try out a variety of initial values for the parameters. Maybe MultiStart is an option.

Maybe your model function is incorrectly implemented - you should check whether .* is always the correct choice against * .

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
stress=stress(2:end);
strain=data(1).strain; %change index to change dogbone if necessary
strain = strain(2:end);
% stress=stress/max(abs(stress));
% strain=strain/max(abs(strain));
p_initial=[0.25 0.55 0.85 1.5];
a0=[1 0 0];
%vincoli su zeta imposti da dong ( zeta tra 0 e 1 )
lb=[0.1 0.1 0.1 0.1];
ub=[inf inf inf inf];
%ottimizzazione non lineare con metodo least square ( like Dong )
%options = optimoptions('lsqnonlin', 'Display', 'iter', 'Algorithm','levenberg-marquardt');
options = optimset('MaxFunEvals',10000,'MaxIter',10000);
[p_opt, resnorm] = lsqnonlin(@(p) error_function(p,strain,stress,a0), p_initial,lb,ub,options);

Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

p_opt

p_opt = 1×4

608.9757 0.1000 1.9532 1.7231

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%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
prod = a0 .* a0';  
I4 = sum(sum(C .* prod)); % 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))) * prod .* 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');

●

function error = error_function(p,strain,stress,a0,k)
   
    sigma_calculated = zeros(length(strain), 1); 
  
    for i = 1:length(strain)
    
    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
    prod = 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)) * prod) .* f';
    sigma_calculated(i)= sigma(1,1);                                           %componente xx
    error(i)= sigma_calculated(i) - stress(i);
     
    end
%il calcolo dell'errore va dentro o fuori dal for?
                           % Somma dei quadrati delle differenze 
%   disp(['Iterazione: ', num2str(k), ' - Errore totale attuale: ', num2str(error)]);
end

[Torsten]

@Nicolò

Is this the correct form of the model function for sigma as a function of strain ?

It's overfitted: the term (zeta-1)^2 and thus zeta can be removed from the model function as a fitting parameter.

syms strain k1 zeta k2 c real
a0 = [1 0 0];
f = [strain 0 0; %deformation gradient for each value of strain
     0 1/sqrt(strain) 0; 
     0 0 1/sqrt(strain)];
C = f .* f'; % Calcola il tensore di Cauchy-Green
product = a0 .* a0';  
I4 = sum(sum(C .* product)); % 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))) * product .* f';
sigma(1,1)

ans = 

[Sam Chak]

Hi @Nicolò

I agreed with @Torsten to double-check the correct form of the model function.

However, according to Zioupos and Barbenel (1994), in their work titled "Mechanics of Native Bovine Pericardium. I. The Multiangular Behaviour of Strength and Stiffness of the Tissue," the graph appears as follows.

If you wish to find a function that describes only the monotonically increasing part, you may consider using non-inverse symmetric sigmoidal (S-shaped) functions, such as the Gompertz distribution or the Soboleva modified hyperbolic tangenthttps://en.wikipedia.org/wiki/Soboleva_modified_hyperbolic_tangent.

You can also stick with the exponential function, but you need to modify it.

load("data_11122.mat")
clearvars -except data
y   = data(1).stress;
y   = y/max(y);
x   = data(1).strain';

Eqn = 'p1*exp(-p2*x)*(p3 + p4*exp(p2*x) + p5*x + p6*x^2 + p7*x^3 + p8*x^4 + p9*x^5)';

fo = fitoptions('Method', 'NonlinearLeastSquares',...
               'Lower', [0, 30, -1, 0, -50, 250, -5e+04, 5e+05, -5e+06], ...
               'Upper', [2, 50,  0, 1, -20, 350, -4e+04, 6e+05, -4e+06], ...
               'StartPoint', [1, 40, -0.5, 0.5, -30, 300, -4.5e+04, 5.5e+05, -4.5e+06]);

ft = fittype(Eqn, 'options', fo);

[yfit, gof] = fit(x, y, ft)

yfit = 
     General model:
     yfit(x) = p1*exp(-p2*x)*(p3 + p4*exp(p2*x) + p5*x + p6*x^2 + p7*x^3 + p8*x^4 
                    + p9*x^5)
     Coefficients (with 95% confidence bounds):
       p1 =       1.744  (-4.279e+04, 4.279e+04)
       p2 =       38.89  (38.56, 39.22)
       p3 =     -0.6848  (-1.681e+04, 1.68e+04)
       p4 =      0.6956  (-1.707e+04, 1.707e+04)
       p5 =         -34  (-8.344e+05, 8.343e+05)
       p6 =       348.3  (-8.546e+06, 8.547e+06)
       p7 =  -4.353e+04  (-1.068e+09, 1.068e+09)
       p8 =   5.619e+05  (-1.379e+10, 1.379e+10)
       p9 =  -4.531e+06  (-1.112e+11, 1.112e+11)

gof = struct with fields:
           sse: 0.0230
       rsquare: 0.9995
           dfe: 416
    adjrsquare: 0.9995
          rmse: 0.0074

plot(yfit, x,  y)
grid on, xlabel('Strain'), ylabel('Stress')
legend('Normalized Data', 'Fitted Function')