Problem with nonlinear curve fitting (lsqcurvefit and Multistart)

Hello all,

Actually, I am dealing with an issue that I can't get desirable results out of my curve-fitting problem. I used lsqcurvefit in the first place. Then, I added MultiStart to my code as well but my result didn't change. I don't know why MATLAB gives the unknown model coefficient the same value as I entered for MATLAB. I bring my code in the following along with the results. I would really appreciate it if you can help me fix this problem and get better results:

clc
clear all;
close all;
%xdata
T = 294:0.1:361;
T=T';
%ydata
for i=1:size(T,1)
f1(i)=  0.2099*exp(-((T(i)-340.2)/2.439).^2);
end
f1=f1';
%Initial value for unknown coefficient
B0 =[-6.1e+04];
%find the local fit using lsqcurvefit
options = optimoptions('lsqcurvefit','FiniteDifferenceType','central','OptimalityTolerance',1e-12,'FunctionTolerance',1e-12, 'MaxFunctionEvaluations',1500);
 lb = [-Inf];
 ub = [Inf];
 [B,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@firstorderDSC1,B0,T,f1,lb,ub,options);
 
% Set up the problem for MultiStart
problem = createOptimProblem('lsqcurvefit','x0',B0,'objective',@firstorderDSC1,...
    'lb',lb,'ub',ub,'xdata',T,'ydata',f1);
ms = MultiStart('PlotFcns',@gsplotbestf);
[xmulti,errormulti] = run(ms,problem,50);
%plot the result
figure;
plot(T,f1,T,firstorderDSC1(xmulti,T),'LineWidth',3)
legend('Data','Fitted result')

Answer 1

[Ameer Hamza]

I am not sure why you wrote the function firstorderDSC1 like that, but reviewing it briefly, it appears that the ode function you have written inside it to fit the data is wrong. Check the differences of your code with this code

clc
clear all;
close all;
%xdata
T = 294:0.1:361;
T=T';
%ydata
for i=1:size(T,1)
    f1(i)=  0.2099*exp(-((T(i)-340.2)/2.439).^2);
end
f1=f1';
%Initial value for unknown coefficient
B0 =rand(1,3);
%find the local fit using lsqcurvefit
options = optimoptions('lsqcurvefit','FiniteDifferenceType','central','OptimalityTolerance',1e-12,'FunctionTolerance',1e-12, 'MaxFunctionEvaluations',1500);
lb = [-Inf];
ub = [Inf];
[B,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@firstorderDSC1,B0,T,f1,lb,ub,options);

% Set up the problem for MultiStart
problem = createOptimProblem('lsqcurvefit','x0',B0,'objective',@firstorderDSC1,...
    'lb',lb,'ub',ub,'xdata',T,'ydata',f1);
ms = MultiStart('PlotFcns',@gsplotbestf);
[xmulti,errormulti] = run(ms,problem,50);
%plot the result
figure;
plot(T,f1,T,firstorderDSC1(xmulti,T),'LineWidth',3)
legend('Data','Fitted result')

function dndt2 = firstorderDSC1(B, T)
%my differential equation
    function denaturate = firstorderDSC(T, n)
        a= (exp(183.8)./1.5); %(1/s)
        dndT=B(1)*exp(-((T-B(2))/B(3)).^2);
        denaturate=dndT;
    end

%solving my differential equation
options = odeset('AbsTol', 1e-8,'RelTol',1e-8);%,'OutputFcn',@odeplot);
[temp,num] = ode23s(@firstorderDSC, T, 0, options);
dndt2 = firstorderDSC(temp,num);
end

Both lines overlap.

[Faezeh Manesh]

Actually, my differential equation is based on my model and it is true. You substituted a Gaussian distribution function with my differential equation which is different. In my equation T is in the denominator.

[Ameer Hamza]

Can you share the mathematical statement of your problem?

[Faezeh Manesh]

Sure. Here is my differential equation.

I have some experimental data (T,dn/dT) and I am trying to fit the result of the above differential equation to my experimental data. Here, instead of experimental data I have (T,f1) which are my xdata and ydata.

[Ameer Hamza]

Is there a reason that you believe that the data generated by the following equation

f1(i)=  0.2099*exp(-((T(i)-340.2)/2.439).^2);

will definitely fit this differential equation?

[Faezeh Manesh]

In the first place, I tried to solve my differential equation and fit the result using gaussian functions and I used curve-fitting tool to do that which is as follows:

So, I am sure that a Gaussian fit to the result of my differential equation. In this code, I am trying to do the opposite. I gave the Gaussian function and I am trying to fit my differential equation to it. Because my ultimate goal is to fit my experimental data using the result of this differential equation.

[Ameer Hamza]

What values of 'a' and 'deltaE/R' you used when creating these curve?

[Faezeh Manesh]

a= (exp(183.8)/1.5);

-deltaE/R= B(1) =-6.2662e+04;

With this value I get the exact same curve as cftool. But I want MATLAB find this value itself.

[Ameer Hamza]

Yes, this makes sense. But this problem is a tough cookie. I have tried several functions to solve this problem using B0=[-6.1e+04], but all failed. Even the functions from the global optimization toolbox didn't work. I guess the problem is caused by the scaling of the terms in the objective function. 'a' is of order 1e80, and the term multiplying with it is of order 1e-80. A tiny change in the parameters can cause a huge difference in the output, which makes it difficult for the optimizer to find an optimal solution. On the contrary, it can also happen that the input parameter changes quite a lot without any effect on the output of ODE, making the optimizer think that it has already reached the optimal solution.

However, I found that if you bring B0 a bit close to the optimal solution, say B0=-6.15e+04, even lsqcurvefit will succeed, even without Multisearch. So, you may need to find a good way to provide an initial guess to solve it with functions available in MATLAB.

[Faezeh Manesh]

Thank you so much for your answer. What are the functions that can provide and initial guess to this problem?

Actually, my ultimate goal is to use multiple of these terms (the solution to these differential equation) to fit to my experimental data. My experimental data are as follows (.csv file is attached) :

I want to fit the above curve by superposing multiple terms came from the differential equation (with differnt a and B). Do you have any idea that how I can do that.

[Ameer Hamza]

Finding a good initial guess can really be considered a theoretical problem. You need to analyze your ODEs to see what type of initial guess will lead to a good solution. For example, you can see that the solution to your ODE has a global peak when you plot T against dn/dT. Is there a way to write down the location of that peak in terms of a and deltaE/R? If yes, you can see that if the peak in the dataset to find the good initial guess for a and deltaE/R (B in your code).

Also, I don't see the reason why you need to use numbers of such large magnitude, e.g., a~1e79. Is there some way to rewrite this system of ODE such that this can be avoided?

In your code, you have fixed the value of parameter 'a'. I have modified your code to tell the optimization algorithm to find both parameters. I have also placed bound on the value of 'a', after some trial and error, I found that it will be useful. Now the solution converges sometime, but the performance is still not reliable enough.

clc
clear all;
close all;
%xdata
T = 294:0.1:361;
T=T';
%ydata
f1 =  0.2099*exp(-((T-340.2)/2.439).^2);
%Initial value for unknown coefficient
B0 =[-6.1e+04 3];
%find the local fit using lsqcurvefit
options = optimoptions('lsqcurvefit','FiniteDifferenceType','central','OptimalityTolerance',1e-12,'FunctionTolerance',1e-12, 'MaxFunctionEvaluations',1500);
lb = [-Inf 1];
ub = [Inf 10];
[B,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@firstorderDSC1,B0,T,f1,lb,ub,options);

figure;
plot(T,f1,T,firstorderDSC1(B,T),'LineWidth',3)
legend('Data','Fitted result')

function dndt2 = firstorderDSC1(B, T)
    %my differential equation
    function denaturate = firstorderDSC(T, n)
        dndT=B(2)*1e79*exp(B(1)./T).*(1-n);
        denaturate=dndT;
    end

    %solving my differential equation
    options = odeset('AbsTol', 1e-8,'RelTol',1e-8);%,'OutputFcn',@odeplot);
    [temp,num] = ode23s(@firstorderDSC, T, 0, options);
    dndt2 = firstorderDSC(temp, num);
end

[Faezeh Manesh]

Thank you so much. You gave me very good ideas.

[Ameer Hamza]

I am glad to be of help.