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
