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)]);
endI 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 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 curve fitting, nonlinear, soft tissue MATLAB Answers — New Questions
