Choosing initial values for fitting experimental data through fminsearch

Hi,

I am trying to find reaction rate constant (k) and order (alpha) for my experimental data of Gas_volume vs. Reaction time.

I used both fminsearch with the help of the general rate law equation. The code is running properly, but the result of the regression does not make sense for our application. I can see that changing the initial guess can have an effect on the final answer, so I wanted to ask what could possibly be the problem? and how to choose the initial guess wisely?

I have checked all the input data multiple times, and there are no issues there. Also, the final values I am expecting should be in the range of [0.0001-0.1 0.99-3]. The value of the first element is compatible with the expectation, but the second value is giving 14, which is way beyond the expected.

Below is the code I am using for the fminsearch:

%Q is the sum of the square residuals 
Q=@(p) sum((C_FA -((((C_FA0)^(1-p(2))-((1-p(2))*p(1)*RxnTime))).^(1/(1-p(2))))).^2);
opt=optimset('fminsearch');
opt.StepTolerance=1E-8;
opt.OptimalityTolerance=1E-14;
opt.FunctionTolerance=1E-14;
%Minimize the sum of the square residuals
P=fminsearch(Q, [0.001 10]);

Your help is appreciated.

Best,

Ayman

[Matt J]

the final values I am expecting should be in the range of [0.0001-0.1 0.99-3].

You cannot allow the search interval for p(2) to contain 1, because the function has a discontinuity there.

[Ayman]

Hi Matt,

Thanks for the reply. This might be the problem behind the strange values for the output.

I am planning to overcome that by adding a nonlinear constraint to the function to have the second value of p not equal to 1. Would that make sense? Or do you suggest a different approach?

Thanks again.

Answer 1

[Star Strider]

It would help to have the data.

If ‘p(2)’ is supposed to be in the range of about 1, start it near there. Nonlinear parameter estimation procedures are sensitive to the initial parameter estimates, and they may not find the values you want if they are initially started farther away.

Perhaps this:

P=fminsearch(Q, [0.001 1]);

or something similar instead.

EDIT — (8 Feb 2024 at 18:16)

I did not initially see the discontinuity that Matt J pointed out. Start with something else close instead, perhaps 0.8 or 2 instead of 1 for the initial ‘p(2)’ estimate.

.

[Ayman]

Thanks for the reply Strider!

I agree that this continuity could be the reason behind the strange data, and as I asked Matt above, I would also appreciate if you have any suggestion on how to overcome that.

Also, attached is the actual data which I am trying to fit using the equation:

Thanks again!

[Star Strider]

Regardless of where I start the second parameter, I get the same result.

Please describe in a bit more detail what you are doing, and specifically the function you want to fit to the data.

T1 = readtable('Dataforfit.xlsx')

T1 = 25×3 table
    RxnTime     C_FA      C_FA0
    _______    _______    _____

        0        0.868    0.868
       15       0.8524      NaN
       30      0.84348      NaN
       45      0.83679      NaN
       60      0.82899      NaN
       75      0.82342      NaN
       90      0.81673      NaN
      105      0.81004      NaN
      120       0.8067      NaN
      135      0.80113      NaN
      150      0.79555      NaN
      165      0.79221      NaN
      180      0.78886      NaN
      195      0.78441      NaN
      210      0.78218      NaN
      225      0.77883      NaN

RxnTime = T1{:,1};
C_FA = T1{:,2};
C_FA0 = T1{1,3};

%Q is the sum of the square residuals 
Q=@(p) sum((C_FA -((((C_FA0)^(1-p(2))-((1-p(2))*p(1)*RxnTime))).^(1/(1-p(2))))).^2);
objfcn = @(p,RxnTime) ((((C_FA0)^(1-p(2))-((1-p(2))*p(1)*RxnTime))).^(1/(1-p(2))));
opt=optimset('fminsearch');
opt.StepTolerance=1E-8;
opt.OptimalityTolerance=1E-14;
opt.FunctionTolerance=1E-14;
%Minimize the sum of the square residuals
P{1} = fminsearch(Q, [0.001 0.01])

P = 1×1 cell array
    {[0.0065 13.8417]}

P{2} = fminsearch(Q, [0.001 1.1])

P = 1×2 cell array
    {[0.0065 13.8417]}    {[0.0065 13.8417]}

figure
plot(RxnTime, C_FA, '.')
hold on
plot(RxnTime, objfcn(P{1},RxnTime), '-r')
hold off
grid

.