I am coding a equation for curve fitting, but I am facing inaccuracy. I am not getting Elastic plateau in the curve. Can anyone guide me how I can address the issue.
113 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
clc; close all;
%% ---- Load data --
file = 'data.xlsx';
T = readtable(file); % expects columns: Strain, Stress (in GPa)
disp(T.Properties.VariableNames)
epsl = T.Strain_mm_mm_(:); % ensure column vector
sig = T.Stress_GPa_(:) * 1000; % convert GPa -> MPa
% keep only finite rows
good = isfinite(epsl) & isfinite(sig);
epsl = epsl(good);
sig = sig(good);
% (optional) trim extreme epsilon to avoid 1/(1-eps) blow-up in model
epsl = min(epsl, 0.981);
%% ---- Known strain rates -------------------------------------------------
edot = 0.000617; % test strain rate [s^-1]
edot0 = edot; % reference strain rate [s^-1]
%% ---- Initial guesses & bounds ------------------------------------------
% [A, E, B, m, n, a, b]
p0 = [0.4, 2000, 0.3, 0.95, 2.0, -0.05, 0.150]; % initial params (MPa-based)
lb = [-0.1, 1.2, -.01, -1, 1, 1, 0.01 ];
ub = [2, 20, 1, 5, 5, 4.0, 0.5 ];
%% ---- Fit with lsqcurvefit (Optimization Toolbox) -----------------------
model = @(p,eps) constitutive_model(eps, edot, p(1), p(2), p(3), p(4), p(5), p(6), p(7), edot0);
try
opts = optimoptions('lsqcurvefit','Display','iter','MaxIterations',1000);
[pfit, ~] = lsqcurvefit(model, p0, epsl, sig, lb, ub, opts);
catch ME
warning('lsqcurvefit unavailable (%s). Falling back to fminsearch (unconstrained).', ME.identifier);
obj = @(p) sum((model(p,epsl) - sig).^2);
pfit = fminsearch(obj, p0);
end
sigma_fit = model(pfit, epsl);
%% ---- Goodness of fit (R^2, adj. R^2, RMSE) -----------------------------
SSE = sum((sig - sigma_fit).^2);
SST = sum((sig - mean(sig)).^2);
R2 = 1 - SSE/SST;
n = numel(sig);
k = numel(pfit); % number of fitted parameters
adjR2 = 1 - (1 - R2) * (n - 1) / max(n - k - 1, 1);
RMSE = sqrt(SSE / n);
% Print to console
fprintf('\nGoodness of fit:\n');
fprintf('R^2 = %.6f\n', R2)
fprintf('Adj R^2 = %.6f\n', adjR2);
fprintf('RMSE = %.6f MPa\n', RMSE)
%% ---- Plot and Save ------------------------------------------------------
fig = figure; hold on;
plot(epsl, sig, 'k-', 'MarkerFaceColor','k','MarkerSize',2, 'DisplayName','Experimental');
plot(epsl, sigma_fit,'r-', 'LineWidth',1, 'DisplayName','Fitted model');
xlabel('Strain \epsilon');
ylabel('Stress \sigma_d (MPa)');
title(['Fitted Constitutive Model at \epsilon-dot = ', num2str(edot), ' s^{-1}']);
legend show; grid on;
% Add a readable stats box
ax = gca;
xL = xlim(ax); yL = ylim(ax);
txt = {sprintf('R^2 = %.4f', R2), ...
sprintf('Adj R^2 = %.4f', adjR2), ...
sprintf('RMSE = %.3g MPa', RMSE)};
text(xL(1) + 0.02*range(xL), yL(2) - 0.04*range(yL), txt, ...
'VerticalAlignment','top', 'BackgroundColor','w', ...
'Margin',6, 'EdgeColor',[0.1 0.1 0.1], 'FontWeight','bold');
% Save the figure
saveas(fig, 'fitted_model_plot.png'); % PNG
% For crisper output you can also use:
% exportgraphics(fig,'fitted_model_plot.png','Resolution',300);
%% ---- Print parameters ---------------------------------------------------
fprintf('\nFitted Parameters:\n');
fprintf('A = %.4f MPa\n', pfit(1));
fprintf('E = %.4f MPa\n', pfit(2));
fprintf('B = %.4f MPa\n', pfit(3));
fprintf('m = %.4f\n', pfit(4));
fprintf('n = %.4f\n', pfit(5));
fprintf('a = %.4f\n', pfit(6));
fprintf('b = %.4f\n', pfit(7));
%% ===================== Local functions (MUST be at file end) =====================
function sigma_d = constitutive_model(epsilon, edot, A, E, B, m, n, a, b, edot0)
% Vector-safe guards
epsilon = min(max(epsilon, 0.001), 0.99); % keep in [0,1)
one_me = max(1 - epsilon, eps);
% Static part
term1 = A .* (1 - exp(-(E./A) .* epsilon .* (1 - epsilon))).^m;
term2 = B .* (epsilon ./ one_me).^n;
sigma_static = term1 + term2;
% Strain-rate factor
rate_factor = 1 + (a + b .* epsilon) .* log(edot ./ edot0);
sigma_d = sigma_static .* rate_factor;
end

0 Commenti
Risposte (1)
Matt J
il 2 Set 2025 alle 0:58
Modificato: Matt J
il 2 Set 2025 alle 1:12
You need a better initial guess. Also, you should exclude the data from the linear elastic region, since clearly that is not captured by your model equation.
file = 'data.xlsx';
[epsl,sig]=readvars(file) ;
sig=sig*1000;% convert GPa -> MPa
% keep only finite and mid-range data
good = isfinite(epsl) & isfinite(sig) & epsl>0.15 & epsl<0.981;
epsl = epsl(good);
sig = sig(good);
%% ---- Known strain rates -------------------------------------------------
edot = 0.000617; % test strain rate [s^-1]
edot0 = edot; % reference strain rate [s^-1]
%% ---- Initial guesses & bounds ------------------------------------------
% [A, E, B, m, n, a, b]
p0 = [0.04, 100, 0.03, 0.95, 2.0, -0.05, 0.150]; % initial params (MPa-based)
lb = [-0.1, 1.2, -.01, -1, 1, 1, 0.01 ];
ub = [2, 20, 1, 5, 5, 4.0, 0.5 ];
%% ---- Fit with lsqcurvefit (Optimization Toolbox) -----------------------
model = @(p,eps) constitutive_model(eps, edot, p(1), p(2), p(3), p(4), p(5), p(6), p(7), edot0);
try
opts = optimoptions('lsqcurvefit','Display','iter','MaxIterations',1000);
[pfit, ~] = lsqcurvefit(model, p0, epsl, sig, lb, ub, opts);
catch ME
warning('lsqcurvefit unavailable (%s). Falling back to fminsearch (unconstrained).', ME.identifier);
obj = @(p) sum((model(p,epsl) - sig).^2);
pfit = fminsearch(obj, p0);
end
sigma_fit = model(pfit, epsl);
postAnalysis(epsl,sig,pfit,sigma_fit,edot)
%% ===================== Local functions (MUST be at file end) =====================
function sigma_d = constitutive_model(epsilon, edot, A, E, B, m, n, a, b, edot0)
% Vector-safe guards
epsilon = min(max(epsilon, 0.001), 0.99); % keep in [0,1)
one_me = max(1 - epsilon, eps);
% Static part
term1 = A .* (1 - exp(-(E./A) .* epsilon .* (1 - epsilon))).^m;
term2 = B .* (epsilon ./ one_me).^n;
sigma_static = term1 + term2;
% Strain-rate factor
rate_factor = 1 + (a + b .* epsilon) .* log(edot ./ edot0);
sigma_d = sigma_static .* rate_factor;
end
function postAnalysis(epsl,sig,pfit, sigma_fit,edot)
%% ---- Goodness of fit (R^2, adj. R^2, RMSE) -----------------------------
SSE = sum((sig - sigma_fit).^2);
SST = sum((sig - mean(sig)).^2);
R2 = 1 - SSE/SST;
n = numel(sig);
k = numel(pfit); % number of fitted parameters
adjR2 = 1 - (1 - R2) * (n - 1) / max(n - k - 1, 1);
RMSE = sqrt(SSE / n);
% Print to console
fprintf('\nGoodness of fit:\n');
fprintf('R^2 = %.6f\n', R2)
fprintf('Adj R^2 = %.6f\n', adjR2);
fprintf('RMSE = %.6f MPa\n', RMSE)
%% ---- Plot and Save ------------------------------------------------------
fig = figure; hold on;
plot(epsl, sig, 'k-', 'MarkerFaceColor','k','MarkerSize',2, 'DisplayName','Experimental');
plot(epsl, sigma_fit,'r-', 'LineWidth',1, 'DisplayName','Fitted model');
xlabel('Strain \epsilon');
ylabel('Stress \sigma_d (MPa)');
title(['Fitted Constitutive Model at \epsilon-dot = ', num2str(edot), ' s^{-1}']);
legend show; grid on;
% Add a readable stats box
ax = gca;
xL = xlim(ax); yL = ylim(ax);
txt = {sprintf('R^2 = %.4f', R2), ...
sprintf('Adj R^2 = %.4f', adjR2), ...
sprintf('RMSE = %.3g MPa', RMSE)};
text(xL(1) + 0.02*range(xL), yL(2) - 0.04*range(yL), txt, ...
'VerticalAlignment','top', 'BackgroundColor','w', ...
'Margin',6, 'EdgeColor',[0.1 0.1 0.1], 'FontWeight','bold');
axis padded
end
14 Commenti
Matt J
il 4 Set 2025 alle 23:21
Modificato: Matt J
il 4 Set 2025 alle 23:27
To demonstrate problems, please get into the habit of running the code here in the online environment, as I have been doing. A screen shot of the results alone doesn't tell us anything, because we cannot see what has changed in the way you are running the code.
By "another strain rate", do you mean a different .xlsx data set?
Matt J
il 4 Set 2025 alle 23:52
Modificato: Matt J
il 5 Set 2025 alle 15:42
One thing that might be relevant to your last question. The known edot and edot0 parameters are not doing anything meaningful in your model, and should probably be removed. All they do is change the fitted values of a and b by a factor of log(edot/edot0). Attached is a version which removes edot,edot0 and also includes a routine to auto-guess the 'c' parameter (so, you no longer have to provide an initial guess or that parameter).
file = 'data_Copy.xlsx';
% [A, E, m, n, a, b]
p0 = [0.1, 500, 0.5, 2.0, 0, 0]; % initial params (MPa-based)
% [A, E, m, n, a, b]
lb = [-0.1, 1.2, -1, 1, -inf, -inf ];
ub = [2, 20, 5, 5, +inf, +inf ];
epsmax=0.85; %Maximum epsilon to fit
doFit(file, p0,lb,ub,epsmax)
%%
file = 'data.xlsx';
% ---- Initial guesses & bounds ------------------------------------------
% [A, E, m, n, a, b]
p0 = [0.9,500, 0.5, 2.0, 0, 0]; % initial params (MPa-based)
% [A, E, m, n, a, b]
lb = [-0.1, 1.2, -1, 1, -inf,-inf ];
ub = [2, 20, 5, 5, +inf, +inf ];
epsmax=0.85; %Maximum epsilon to fit
doFit(file, p0,lb,ub,epsmax)
Vedere anche
Categorie
Scopri di più su Simulink Report Generator in Help Center e File Exchange
Prodotti
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!