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

3 views (last 30 days)
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')
  8 Comments
Torsten
Torsten on 5 Apr 2019
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).

Sign in to comment.

Accepted Answer

Alex Sha
Alex Sha on 4 Jun 2019
Edited: Alex Sha on 4 Jun 2019
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
c171.jpg
c172.jpg
c173.jpg

More Answers (2)

Alex Sha
Alex Sha on 28 Apr 2019
Please upload your data,I may try.

John D'Errico
John D'Errico on 5 May 2019
Edited: John D'Errico on 5 May 2019
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-')
untitled.jpg
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.
  2 Comments
Kaouthar Kaouthar
Kaouthar Kaouthar on 31 May 2019
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
Kaouthar Kaouthar on 31 May 2019
Hello, I did not understand your proposed code already. can you help me understand how you defined sp1, sp2 and sp3 please?
Thank you

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!