Parameter estimation with lsqcurvefit and ODE45

Ok, I am must be doing sometime wrong and it is time to get help. I am trying to estimate solid state kinetic parameters for a set of data using the Avrami-Erfoveev which is f(a) = n(1 - a)(-ln(1-a))^(n-1)/n where a (or alpha) is the fractional conversion and n is the reaction order. The rate law for this reaction is then da/dt = A*exp(Ea/R*T)*f(a). I am trying to estimate the A, Ea and n parameters for three set of differential equations.

I have created the following matlab nested functions largely from answers provided on members.

% data = importfile('sample.mat'); %imports data from a data file
% alphas = [data.ah data.ac data.al data.a];
%
% time = data.t; Time (sec)
% T = data.T; Temperature (K)
%
%   
% function C=Avramikinetics(beta,time)
%     a0 = [alphas(1,1); alphas(1, 2); alphas(1,3); alphas(1,4)]; % Initial conditions
%
%    R = 8.314; gas constant
%
%    ode_options = odeset('RelTol',1e-8, 'AbsTol', 1e-08, 'Stats','on'); % Changed per Torsten's suggestion
%    [T,Cv]=ode15s(@DifEq,time,a0, ode_options); % Changed per Torsten's suggestion
%
%    function da=DifEq(time,alphas)
%      dadt=zeros(4,1);
%
%      dadt(1) = beta(1,1) .* exp(-beta(1,2) ./ (R *T(1))) .* beta(1,3) * (1- alphas(1)) .* (-log(1- alphas(1))) .^ ((beta(1,3) - 1)/beta(1,3));
%      dadt(2) = beta(2,1) .* exp(-beta(2,2) ./ (R *T(1))) .* beta(2,3) * (1- alphas(2)) .* (-log(1- alphas(2))) .^ ((beta(2,3) - 1)/beta(2,3));
%      dadt(3) = beta(3,1) .* exp(-beta(3,2) ./ (R *T(1))) .* beta(3,3) * (1- alphas(3)) .* (-log(1- alphas(3))) .^ ((beta(3,3) - 1)/beta(3,3));
%      dadt(4) = dadt(1) + dadt(2) + dadt(3);      
%      da=dadt;
%     end
% C=Cv;
% end
%
% beta0= [1e+30 2e+05 1; 1e+30 2e+05 1;5e+13 1e+05 1]; %initial guesses
% lower = [0, 0, 0; 1e+29, 0, 0; 1e+11, 0, 0]; %lower boundary for paramaters
% upper = [1e+035, 1e+06, 3; 1e+035, 1e+06, 3; 1e+018, 1e+06, 3]; %upper boundary for paramaters
%
% options = optimoptions('lsqcurvefit','Algorithm','trust-region-reflective');
% options.OptimalityTolerance = 1.0000e-016;
% options.FiniteDifferenceStepSize = 1; % Changed per Alan Weiss's suggestion. The default was 1e-03 I believe
% options.FiniteDifferenceType = 'central';
% options.Diagnostics = 'on';
% options.Display = 'iter';
% options.FunctionTolerance = 1.0000e-016;
% options.FunValCheck = 'on';
% options.MaxIterations = 1000;
%options.StepTolerance = 1.0000e-016;
%
% [beta, resnorm, residual, exitflag, output, lambda, jacobian] = lsqcurvefit(@Avramikinetics,beta0,time,alphas, lower, upper);

When I run this code however, it never converges. Too often, it spits out the unchanged initial estimates as final. I will greatly appreciate any help. I have attached a sample dataset.

Thanks

[Torsten]

1. The Arrhenius term must be coded as exp(-Ea/(R*T)), not as exp(-Ea/R*T).

2. The dimensions of the arrays lower and upper must be the same as those of beta0, thus 3x3.

3. Since the equations for the three a-curves are independent, I'd first try to fit only one of them.

4. In a first step, I'd calculate the a-curve for your initial guess of the parameters and compare with your data curve to see if they are quite close to each other.

Best wishes

Torsten.

[Ahn]

Torsten,

Thank you for the comments and suggestions which I have implemented with few more addition to make sense of what lsqcurvefit is doing behind the scene. Regarding your comment#3: While the first three equations are independent of each other, all three are tied to the fourth equations as shown in modified code. The overall progression of the reaction (reaction 4) is a function of three parallel reactions (reactions 1-3).

Regarding your comment#4: I have improved on the initial guesses by using the Avramikinetics function and initial guesses to get a sense of what reasonable numbers to provide as initial guesses.

However, despite these changes, I am still running into the same challenges. In particular, lsqcurvefit seems to only go through two iterations before it spits back out the same initial guesses provided. Am I not using lqscurvefit properly or is my implementation simply off? Btw, I have based my general approach based off a post here by Star Strider.

Thanks

[Torsten]

I suggest you switch to ODE15S as ODE solver and use a smaller tolerance than 1e-5 (e.g. 1e-8 for both RelTol and AbsTol).

Best wíshes

Torsten.

[Ahn]

I changed to ODE15s reduced the tolerance. Additionally, I changed FiniteDifferenceType setting from 'forward' to 'central'. Finally, I changed the larger finite differences from 1e-03 to 1 based Alan's suggestion. After those changes, the initial and final estimates are identical. However, unlike in the past, lsqcurvefit iterated about 20. I monitored this by turning Diagnostics on. I am a bit baffled here. I am attaching two plots: #1 was generated using my initial guesses without any lsqcurvefit regression and #2 was generated after lsqcurvefit using those same initial guesses.

[Torsten]

The curves seem to be identical. In the one case, you plot "a", in the other, you plot "1-a".

Best wishes

Torsten.

[Ahn]

I have edited the curves for more clarity. As you can see in #2.png, the simulation curves are significantly off.

[Torsten]

Sorry, I still don't understand completely what your graphics show.

As far as I understand, #1.png shows your experimental curves together with the results of the integration of your four ODEs with an initial guess of the parameters (without using lsqcurvefit).

But what does #2.png show (especially because simulated curves seem to be missing in the graphics) ?

Aside from this:

You will have to go one step back. The usual way to proceed now is as follows:

1. Choose an arbitrary set of parameters.

2. Integrate your ODEs using these parameters.

3. Perturb the results obtained by small quantities.

4. Take these perturbed simulation curves as "experimental" curves and try to recover your parameters now with the help of "lsqcurvefit".

You should take initial guesses for the parameters not identically to those you chose, but not too far away.

This procedure usually shows whether your code works out correctly.

Best wishes

Torsten.

[Ahn]

Torsten,

In #2.png, the simulated curves are one on top of another, at the baseline giving the appearance that they are not present. Effectively, in one iteration, lqscurvefit seems to change the estimates to values that deviate significantly from the initial guesses. Then those estimates are not changed. I wonder what objective function it is minimizing since those estimates would most certainly not minimize the sum of squares.

I will give your latest suggestion a try...

Answer 1

[Alan Weiss]

You might be getting bit by too-small finite differences, which sometimes happens when your objective function is given by the solution to an ODE. The link has solutions you can try, meaning how to set larger finite differences.

Alan Weiss

MATLAB mathematical toolbox documentation

[Ahn]

Alan,

Thanks for the suggestion. Upon implementing your suggestion, the issue persists. See my reply to Torsten for additional details.