Line search algorithm: nonlinear least squares method for fitting experimental data to a model

Hello, I have tried to fit my experimental data with a nonlinear model. But the error between the obtained curves is very high. I tried to solve it with lsqnonlin and lsqcurvefit but the exitflag that I received is 1, so i can not get better results with this model. I find then that in some articles, they used line search method. I am not familiar with these types of codes. Can you, please, help me adapting my code to line search algorithm? Otherwise, any ideas to get better fit.

Thank you in advance.

clc
close all
clear all
figure;
A=importdata('1.xlsx');
xdata = ...
 [A(:,1) A(:,1) A(:,1)];
ydata = ...
 [A(:,2) (A(:,3)) (A(:,4))];
%%
% Create the  model.
fun = @(x,xdata)[(x(1)*exp(-(x(2)*(357.15-x(4)))/(x(3)+357.15-x(4))))...
        ./(1.+(xdata(:,1).*(x(1)*exp(-(x(2)*(357.15-x(4)))./(x(3)+357.15-x(4))))./x(6)).^(1-x(5))),...
        (x(1)*exp(-(x(2)*(367.15-x(4)))./(x(3)+367.15-x(4))))...
        ./(1.+(xdata(:,1).*(x(1)*exp(-(x(2)*(367.15-x(4)))./(x(3)+367.15-x(4))))./x(6)).^(1-x(5))),...
        (x(1)*exp(-(x(2)*(377.15-x(4)))./(x(3)+377.15-x(4))))...
        ./(1.+(xdata(:,1).*(x(1)*exp(-(x(2)*(377.15-x(4)))./(x(3)+377.15-x(4))))./x(6)).^(1-x(5)))];
%%
% Fit the model using the starting point 
% x0 = [1e5;10;10;340;0.2;500];
x0 = [1e4;2;2;300;0.3;6];
options = optimoptions('lsqcurvefit','Algorithm','levenberg-marquardt','Display','iter','TolX', 1e-6, 'TolFun', 1e-6, 'MaxFunEvals', 40000, 'MaxIter', 400);
lb = [1 1 1 1 0 1];
ub = [1e17 100 100 1000 1 1e9];
x =  lsqcurvefit(fun,x0,[xdata(:,1) xdata(:,1) xdata(:,1)],[ydata(:,1) ydata(:,2) ydata(:,3)],lb,ub,options)
%%
% Plot the data and the fitted curve.
loglog(xdata(:,1),fun(x,xdata(:,1)))
 hold on
loglog(xdata(:,1),ydata(:,1),'ko')
hold on
loglog(xdata(:,1),ydata(:,2),'ro')
hold on
loglog(xdata(:,1),ydata(:,3),'go')
legend('Fitted1','Fitted2','Fitted3','Data1','Data2','Data3')
title('Data and Fitted Curve')

[Torsten]

I guarantee that if lsqcurvefit does not succeed fitting your model, you won't succeed with a self-written line search code, either.

[Kaouthar Kaouthar]

Thank you for your reply Torsten. Can you please have a look at this attached article? They used this method for exactly the same problem. Is that correct what I understood?

[Torsten]

Please check the three curves you get if you use the initial parameter values. Are they far away from your experimental results ?

[Kaouthar Kaouthar]

Yes, they are far away the experimental data and they the shape of a line (not a curve similar to the experimental data). I don t know even how to choose the start point. Do you have idea about this?

Thank you again for your quick reply.

[Torsten]

Can't you choose the values from the article ?

[Kaouthar Kaouthar]

I ve tried them but we don t have the same polymer materials. That s why they do not fit well my model.

[Torsten]

Then you should first play around with the parameters.

Without good initial guesses, your fitting procedure won't be successful (independent of the method you try).

[Kaouthar Kaouthar]

Okay, i will try to find a good start point. Thank you for your interest and follow-up.

Answer 1

[Alex Sha]

Hi, Kaouthar, the code changed to fellow, it ensures all parameters are postive, and the label of x_axis is your true xdata. However one thing should be noticed: the parameters are not unique for x1, x2 x3 and x4.

VarConstant c=[357.15,367.15,377.15];
SharedModel;
Parameter x(6)>0;
XAxis = xdata;
Variable xdata,y;
Function  y=(x1*exp(-(x2*(c-x4))/(x3+c-x4)))/(1+(xdata*(x1*exp(-(x2*(c-x4))/(x3+c-x4)))/x6)^(1-x5));
Data;
0.5,0.659,0.869,1.15,1.51,1.99,2.63,3.47,4.57,6.03,7.95,10.5,13.8,18.2,24.1,31.7,41.8,55.2,72.7,95.9,127,167,220,290,383,505,665,878,1160,1530,2010,2650,3500;
24.434,24.336,22.598,19.534,16.388,13.596,11.234,9.282,7.692,6.514,5.649,4.839,4.161,3.602,3.159,2.792,2.499,2.262,2.070,1.916,1.792,1.687,1.601,1.537,1.501,1.473,1.447,1.418,1.387,1.345,1.297,1.228,1.072;
Data;
0.5,0.659,0.869,1.15,1.51,1.99,2.63,3.47,4.57,6.03,7.95,10.5,13.8,18.2,24.1,31.7,41.8,55.2,72.7,95.9,127,167,220,290,383,505,665;
20.583,19.918,18.385,15.886,13.327,11.047,9.116,7.537,6.299,5.326,4.484,3.814,3.249,2.796,2.438,2.140,1.905,1.721,1.567,1.442,1.341,1.259,1.194,1.145,1.115,1.096,1.078;
Data;
0.5,0.659,0.869,1.15,1.51,1.99,2.63,3.47,4.57,6.03,7.95,10.5,13.8,18.2,24.1,31.7,41.8,55.2,72.7,95.9,127,167,220,290,383,505,665;
17.141,16.598,15.950,14.201,12.054,10.024,8.269,6.820,5.642,4.709,3.946,3.326,2.826,2.417,2.088,1.823,1.609,1.438,1.297,1.185,1.094,1.021,0.962,0.919,0.889,0.868,0.851;

Result 1:

Root of Mean Square Error (RMSE): 0.750624723639767

Sum of Squared Residual: 49.0190603893171

Correlation Coef. (R): 0.994414241612346

R-Square: 0.988859683921457

Adjusted R-Square: 0.988000769669663

Determination Coef. (DC): 0.986551039923055

F-Statistic: 342.846071379397

Parameter Best Estimate

-------------------- -------------

x1 50.7166624975047

x2 1.74233039846288

x3 61.0256533126286

x4 370.012581419946

x5 0.34018118177961

x6 14.1379270601172

Result 2:

Root of Mean Square Error (RMSE): 0.750624723639768

Sum of Squared Residual: 49.0190603893173

Correlation Coef. (R): 0.994414241654891

R-Square: 0.988859684006072

Adjusted R-Square: 0.988000769760861

Determination Coef. (DC): 0.986551039936596

F-Statistic: 342.846071879207

Parameter Best Estimate

-------------------- -------------

x1 67.2294239650589

x2 2.02418733773269

x3 52.5281917961558

x4 361.515106529892

x5 0.340181179660682

x6 14.1379272335838

Result 3:

Root of Mean Square Error (RMSE): 0.750624723639767

Sum of Squared Residual: 49.019060389317

Correlation Coef. (R): 0.994414241636979

R-Square: 0.988859683970447

Adjusted R-Square: 0.988000769722583

Determination Coef. (DC): 0.986551039933522

F-Statistic: 342.846071436732

Parameter Best Estimate

-------------------- -------------

x1 743.828260682761

x2 4.42788617741519

x3 24.0130073020205

x4 332.999932457755

x5 0.34018118195002

x6 14.1379269290681

[Kaouthar Kaouthar]

Thank you very much

[kelly]

@Alex Sha I don't understand what you did to resolve the question of @Kaouthar Kaouthar, can you please explain it to me? Thank you.

[anahita]

hello can you just find error of this code. i cant have a plot like the plot you have clc; clear; close all;

data = readmatrix('analysis_1.xlsx', 'Sheet', 'Main');

S = data(:, 6);

mu = data(:, 2);

sigma = data(:, 3);

gamma = data(:, 4);

Y = data(:, 5);

ZR = data(:, 7);

%% Ryan-Joiner mothode

sorted_S = sort(S);

n = length(S);

z_scores = norminv(((1:n) - 0.5) / n, 0, 1);

S_mean = mean(S);

S_std = std(S);

e_i = (sorted_S - S_mean) / S_std;

Rp = sum(e_i .* z_scores) / sqrt(S_std^2 * (n - 1) * sum(z_scores.^2));

disp(['Ryan-Joiner Rp = ', num2str(Rp)]);

%% (Nonlinear Least Squares)

% S = a * mu^b * sigma^c * gamma^d * Y^e * ZR^f

modelFun = @(params, X) params(1) * (X(:,1) .^ params(2)) .* ...

(X(:,2) .^ params(3)) .* (X(:,3) .^ params(4)) .* ...

(X(:,4) .^ params(5)) .* (X(:,5) .^ params(6));

X = [mu, sigma, gamma, Y, ZR];

initial_guess = [1, 1, 1, 1, 1, 1];

params_opt = lsqcurvefit(modelFun, initial_guess, X, S);

a = params_opt(1);

b_coeff = params_opt(2:end);

disp(['a = ', num2str(a)]);

disp(['b = ', num2str(b_coeff(1))]);

disp(['c = ', num2str(b_coeff(2))]);

disp(['d = ', num2str(b_coeff(3))]);

disp(['e = ', num2str(b_coeff(4))]);

disp(['f = ', num2str(b_coeff(5))]);

%% R²

S_Estimate = modelFun(params_opt, X);

SS_res = sum((S - S_Estimate).^2);

SS_tot = sum((S - mean(S)).^2);

R2 = 1 - (SS_res / SS_tot);

disp(['R2 = ', num2str(R2)]);

results_table = table(S, S_Estimate, 'VariableNames', {'Actual_Storage', 'Estimated_Storage'});

coeff_table = array2table([a; b_coeff'], 'VariableNames', {'Value'}, ...

'RowNames', {'a', 'b', 'c', 'd', 'e', 'f'});

r2_table = table(R2, 'VariableNames', {'R2_Value'});

filename = 'matlab_nonlinear_regression1.xlsx';

writetable(results_table, filename, 'Sheet', 'Storage Comparison');

writetable(coeff_table, filename, 'Sheet', 'Regression Coefficients', 'WriteRowNames', true);

writetable(r2_table, filename, 'Sheet', 'R2 Value');

disp('Excel report generated successfully!');

figure;

scatter(S, S_Estimate, 'b', 'filled');

hold on;

plot([min(S), max(S)], [min(S), max(S)], 'r--', 'LineWidth', 2);

hold off;

xlabel('Actual Storage (S)');

ylabel('Estimated Storage (S_{Estimate})');

title(['Nonlinear Regression (R^2 = ', num2str(R2), ')']);

legend('Data Points', '1:1 Line', 'Location', 'Best');

grid on;

Answer 2

[Alex Sha]

Please upload your data，I may try.

[Kaouthar Kaouthar]

i have already insertd my data "1.xlsx"

Answer 3

[John D'Errico]

A bad idea: to write your own code to implement a numerical method, because you cannot make professionally written code to do the same. Writing your own line search code is exactly that - a really bad idea.

Instead, learn why it is that your model fit has problems.

xdata = [A(:,1);A(:,1);A(:,1)];
ydata = [A(:,2);A(:,3);A(:,4)];
fun = @(x,xdata) (x(1)*exp(-(x(2)*(357.15-x(4)))/(x(3)+357.15-x(4))))...
        ./(1.+(xdata(:,1).*(x(1)*exp(-(x(2)*(357.15-x(4)))./(x(3)+357.15-x(4))))./x(6)).^(1-x(5)));
    
lb = [100 1 1 1 0 1];
ub = [5e5 100 100 1000 1 1000];

x0 = [1e4;2;2;300;0.3;6];
options = optimoptions('lsqcurvefit','Algorithm','levenberg-marquardt','Display','iter','TolX', 1e-6, 'TolFun', 1e-6, 'MaxFunEvals', 40000, 'MaxIter', 400);
lb = [1 1 1 1 0 1];
ub = [1e17 100 100 1000 1 1e9];
x =  lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options)

x
x =
          9651.15786422703
          5.32989985066721
          1.96617535975919
          289.839867479332
         0.325125105721047
          15.7701020874321
          

Now, what does the data plot suggest?

loglog(xdata,ydata,'ro',xdata,fun(x,xdata),'b*')
grid on

Ok, we see that the curve does not fit well at the bottom end. This stems from two reasons. your fit was done in terms of a sum of squares, but the plot is show on loglog axes. Down at the bottom end, errors in those points are virtually irrelevant when you compute a straight sum of squares. They are tiny numbers compared to the rest, and the algorithm does not "care" since the errors down there are tiny.

As bad, you are trying to fit all three curves with the same set of parameters. They are very clearly different curves. Wanting to use the same parameters for all three curves seems a bit silly, because the variation between curves is quite large.

I'd suggest that either your model simply does not fit that data, or your starting values for the model are terribly poor. But using a different search method is not the way to fix this problem.

Essentially, the question I would ask is if this model actually means something to you in context? Do those parameters make physical sense? Did you derive this model from first physical principles, or did you just pick some terms that might give the curve some flexibility? Lets see, maybe I'll add an exponent here? So, maybe you started with something that makes some vague physical sense in context, but then you started adding terms and parameters?

This is all important, because what you intend to do with that model will matter. Because each of those three curves are fundamentally different in shape from each other, even while they all have distinct similarities too. So the question strongly becomes, why are you trying to fit all three of those curves with one set of parameters, when your model does not even fit any of them well? Thus a simple set of interpolating splines will do very nicely? Or a simple cubic least squares spline could be used to fit each curve, and nail down the behavior of each of those curves.

For example:

xdata = A(:,1);
ydata = [A(:,2) (A(:,3)) (A(:,4))];

loglog(xdata,ydata,'ro')
hold on
grid on
loglog(xdata,exp(fnval(sp1,log(xdata))),'b-')
loglog(xdata,exp(fnval(sp2,log(xdata))),'g-')
loglog(xdata,exp(fnval(sp3,log(xdata))),'k-')

Now, if the various parameters in your model actually have some meaning in context about the shape of that curve, you could now try to extracct them from the shape of those splines too. And you would have far more success in the process.

[Kaouthar Kaouthar]

Hello John,

I just noticed that you have modified your answer, so sorry for the late to reply and thank you again for your interest.

So, to reply to your questions, please consult these two pdf documents (attached) that treat exactly the same problem as me. Thus, you will undersand why this model? physics behind? why three simultaneously? etc

Thank you again.

[Kaouthar Kaouthar]

Hello, I did not understand your proposed code already. can you help me understand how you defined sp1, sp2 and sp3 please?

Thank you