how to control fitting in 'lsqcurvefit'?

1 回表示 (過去 30 日間)
Monirul Hasan
Monirul Hasan 2018 年 12 月 7 日
回答済み: John D'Errico 2018 年 12 月 7 日
Hi, I have a experimental data fitting code given below. In this code from the differential equation (line 30), the 'lsqcurvefit' function (line 12) giving B(1) and B(2) estimation after fitting. How can like give a condition like B(1)+B(2)=9.6*10^-1 before fitting so based on that condition lsqcurvefit will find the best fit.
Please feel free to check the attachement for the experimental data.
clear all, close all, clc
HighVacData = xlsread('RawData','Sheet1','A2:B1603');% loading raw data from excel file
Timeall = HighVacData(:,1);
Sdataall = HighVacData(:,2);
Time = 0:1:11244.3353200000;
TimeC = Time.';
%parameter estimation (units per ns)
B0 = [2.96e-02,1.45e-3]; %sequence [k_rs]
lb = [0,0];% lower bound estimation
options = optimset('Algorithm', 'trust-region-reflective');
[B,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@ComGenUQFit,B0,Timeall,Sdataall,lb,[],options);%sending data for fitting
F = ComGenUQFit(B,TimeC);% fitted data saving
figure(1)
plot(Timeall,Sdataall,'bo')%raw data
hold on
plot(TimeC,F,'r','linewidth',2)%fit data
grid on
grid minor
xlabel('Time [ns]');
ylabel('Singlet');
legend('Experimental Data','Fit');
set(gca,'xscale','log')
set(gca,'yscale','log')
function S = ComGenUQFit(B,TimeC)% function for fitting
IC_ns = 11000;
function dotn = Mat (TimeC,n)
dotn(1) = -B(1)*n(1)-B(2)*n(1)^2;% equation for singlet Population
end
[TimeC,dotnsolve] = ode45(@Mat,TimeC,IC);% calling ODE solver
S = dotnsolve(:,1);
end

採用された回答

Rik
Rik 2018 年 12 月 7 日
Since you can calculate B(2) based on B(1), your function has actually only a single dependent, so you can adapt your function like this to estimate B(1), and then after fitting, calculate B(2).
function S = ComGenUQFit(B,TimeC)% function for fitting
IC_ns = 11000;
B=[B 0.96-B];%fix B(1)+B(2)=9.6*10^-1
function dotn = Mat (TimeC,n)
dotn(1) = -B(1)*n(1)-B(2)*n(1)^2;% equation for singlet Population
end
[TimeC,dotnsolve] = ode45(@Mat,TimeC,IC);% calling ODE solver
S = dotnsolve(:,1);
end
(of course you should adapt your B0 initial guess as well)

その他の回答 (1 件)

John D'Errico
John D'Errico 2018 年 12 月 7 日
You cannot do so with lsqnonlin or lsqcurvefit, because those tools do not allow general linear equality constraints. Sorry, but not an option. Similarly, nlinfit or the curvefitting toolbox will not allow that class of constraint.
You have several choices. A simple one if to convert the problem to one that fmincon can handle, since it allows equality constraints. That means you need to artifically form andreturn the sum of squares of the residuals as the objective. It will also slightly reduce the final accuracy you could obtain, because tools like lsqnonlin/lsqcurvefit do not need to explicitly square the residuals in the linear algebra. This is a subtle point, but I am sure someone would trip over it one day.
The second option is to use a transformation, reducing the number of parameters in your model.
So in your model, just replace B(2) with the term (9.6*10^-1 - B(1)), in every place B(2) was used. Make sure you remember to renumber the other parameters too, in case B(3) was in the model, etc.
Yes, I expect this will not make you happy, that you don't want to think about changing your model form. People never seem to like that idea.
Finally, with only some little effort, one could write an overlay to tools like lsqnonlin or lsqcurvefit, that would do the model reduction for that class of linear constraint explicitly, then internally call lsqnonlin/lsqcurvefit as required with the reduced model, but returning the result in the original form for the variables. Personally, it seems like more work than it is worth, since it is too easy to do with pencil and paper to just change the model. Yes, I could write the code, or they could have implemented it in lsqnonlin.

製品

Community Treasure Hunt

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

Start Hunting!

Translated by