Initialization values are no changing in an ODE data fitting with lsqnonlin. Is there a way to check initialization values faster? this method takes too long

Good evening everyone,

I am trying to estimate kinetic parameters by fitting experimental data to 4 ODE. I have 8 parameters to estimate and 4 ODEs to resolve. For this i am using ode45 to solve the ODE and lsqnonlin to optimize the kinetic parameters in order to fit experimental data. Seems there is not a syntax problem on my code since it runs properly but the result of the parameters is exactly the same to the initialized values for some of them, even if i change them extremely (from an initialization value of 100 to 1e20) so maybe something inside my code is wrong but since it runs i haven't been able to find it. After optimized, the objective function always gives the same result with no change in variables.

here is the code i am using..

clc;
clear all;

%%Initialization of the kinetic parameters

b0=[10000;10000;10000;10000;...
    100000;100000;100000;100000];

%%Optimization parameters

LB=[1e5,1e5,1e5,1e5,1e5,1e5,1e5,1e5];%Lower bounds of every parameter

UB=[Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf];%Upper bounds of every parameter

options=optimset('Disp','iter',... %Display the iterations
'MaxIter',1e9,...%Maximum number of iterations allowed.
'TolX',1e-20,...%Termination tolerance on the variable
'TolFun',1e-20,...%Termination tolerance on the function value. 
'LargeScale','on',...%Use large-scale algorithm 
'MaxFunEvals',1E10,...%Maximum number of function evaluations allowed
'Algorithm','levenberg-marquardt');

%%Solvers

%lsqnonlin Solve nonlinear least-squares (nonlinear data-fitting) problems

[b,resnorm,residual,exitflag,output,lambda]=lsqnonlin('run',b0,LB,UB,options);

----

function SSE = run(b)

global Acel Aglc1 Aglc2 Ahmf Ecel Eglc1 Eglc2 Ehmf mcel mglc1 mglc2 mhmf 

rslts =[1.7268  1.8862  0.4731  0.1477  1.2394  0.0000;
        1.6150  1.4807  0.7366  0.2122  1.039  0.0000;
        1.6418  1.4119  0.7545  0.2383  1.5249  0.0000;
        1.6472  1.3797  0.7239  0.2594  1.3542  0.0000;
        1.5349  1.2926  0.7018  0.3127  1.2197  0.0000;
        1.7451  1.4772  0.7824  0.4353  1.2928  0.3001;
        1.7338  1.4901  0.1134  0.4925  1.2382  0.3252;
        1.1987  1.0593  0.5487  0.4322  0.8575  0.2760;
        1.3623  1.2195  0.6256  0.5848  0.9279  0.3717;
        1.3937  1.3277  0.2874  0.6387  0.9099  0.4181;
        1.4624  1.4136  0.3416  0.7415  0.9620  0.4752;
        1.4282  1.3930  0.5202  0.7877  0.9321  0.4830;
        0.9200  0.9212  0.3902  0.5670  0.6456  0.3580;
        1.2116  1.2165  0.5641  0.8967  0.9217  0.5026;
        1.0532  1.1171  0.5621  0.9196  0.8419  0.4828;
        0.9515  1.1064  0.6549  1.0482  0.8897  0.5479;
        0.6938  0.9480  0.6690  1.0692  0.8134  0.5126];

concGlc = rslts(:,1);

concLA  = rslts(:,4);

concHMF = rslts(:,5);

%%parameters to be found

Acel  = b(1);

Aglc1 = b(2);

Aglc2 = b(3);

Ahmf  = b(4);

Ecel  = b(5);

Eglc1 = b(6);

Eglc2 = b(7);

Ehmf  = b(8);

      function f = odefun(t,x)

      CEL = x(1);
      GLC = x(2);
      HMF = x(3);

        %%Constants

         R=8.314;         % Gas constant
         Temp=175+273;        % Experiemntal Temperature (K)
         Acd= 0.00533;        % Expermiental acid concentration(%w/w) 

        %%kinetic model

         kcel=Acel*exp(-Ecel/(R*Temp))*Acd;%^mcel;
         kglc1=Aglc1*exp(-Eglc1/(R*Temp))*Acd;%^mglc1;
        kglc2=Aglc2*exp(-Eglc2/(R*Temp))*Acd;%^mglc2;
        khmf=Ahmf*exp(-Ehmf/(R*Temp))*Acd;%^mhmf;

         rcel=kcel*CEL;

        rglc1=kglc1*GLC;

        rglc2=kglc2*GLC;

        rhmf=khmf*HMF;

        %%Differential equations

        dCELdt = -rcel;

        dGLCdt = rcel-rglc1-rglc2;

        dHMFdt = rglc1-rhmf;

        dLAdt = rhmf;

        f = [dCELdt,dGLCdt,dHMFdt,dLAdt]';
        end

tspan = [1 2 3 4 5 7 9 11 13 15 17 21 25 30 40 50 60];

y0 =[0.8029 1.7268 1.2394 0.1477];

[v y] = ode45(@odefun , tspan, y0);

ymeas= [concGlc concHMF concLA];

yest = [y(:,2) y(:,3) y(:,4)];

   err=yest-ymeas;       

   SSE = abs(err);     % for lsqnonlin
  end

I would certainly appreciate all your help. If you find an error in my code or if you know a way to make it better I would appreciate if you say so.

Thanks in advvance!

[Nath]

Hi There

I got similar problem any suggestion on how to solve it. Did you get any better answer?

Thankx

Answer 1

[Alan Weiss]

Without reading your code in detail, I can suggest the following:

• Don't set tolerances below eps ~ 2e-16. See the documentation, where a minimum of 1e-14 is suggested.
• Don't give the UB vector, since it is all Inf anyway. Set UB to [].
• Don't set the levenberg-marquardt algorithm, because that algorithm does not handle bound constraints.
• Consider scaling your problem; your lower bound values are pretty large.
• Whether or not you scale, set (much) larger finite differences for gradient estimation, as explained in the documentation.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation

[nicolin2001]

Ok thank you. I just realized the LB problem. However, what do you think a larger finite differences set would be appropriate?

I passed prom 1e-3 as default to 10 but i'm not sure if would it be. i still need accuracy and definitely i am not used to optimization issues.

Thanks once again