Final Point Equals Initial Point when trying to optimize

The initial guesses for A0 are equal to what has been calculated for A. How can i trouble shoot?

I tried changing the initial values, changing the tolerance, changing functions... non seems to be working.

[Torsten]

Translating your code, you use as time interval for integration [299:1307.99] and insert the actual time T during integration in the expressions exp(-Ei/(R*T)). This means that the temperature of your system rises from 299 to 1307.99 K linearly with time.

Is this really what you intend to do ?

And don't name your array of experimental values "exp". It can easily come into conflict with the MATLAB exponential function exp(x) for e^x.

[Rosalind]

Yes, temperature linearly increases. I changed the name from exp to experiment but still have the same problem.

Answer 1

[Torsten]

Here is a code where A does not equal A0:

A0=1e-3*[1.5*(10^7),1*(10^6),1*(10^6)]

A0 = 1×3

15000 1000 1000

sum(fit(A0).^2)

ans = 0.0945

options = optimset('TolFun',1e-7);
[A,resnorm,residual,exitflag,output,lambda,Jacobian] = lsqnonlin(@fit,A0)

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

A = 1×3

1.0e+03 * 2.2446 1.3592 0.3087

resnorm = 0.0404

residual = 4×1

0.1025 0.0994 0.0999 0.1001

exitflag = 1

output = struct with fields:
    firstorderopt: 2.9048e-07
       iterations: 17
        funcCount: 72
     cgiterations: 0
        algorithm: 'trust-region-reflective'
         stepsize: 393.3516
          message: 'Local minimum found.↵↵Optimization completed because the size of the gradient is less than↵the value of the optimality tolerance.↵↵<stopping criteria details>↵↵Optimization completed: The first-order optimality measure, 2.904751e-07,↵is less than options.OptimalityTolerance = 1.000000e-06.'

lambda = struct with fields:
    lower: [3×1 double]
    upper: [3×1 double]

Jacobian = 
   1.0e-03 *

   (1,1)      -0.0235
   (2,1)       0.0452
   (3,1)      -0.0126
   (4,1)      -0.0091
   (1,2)      -0.0467
   (2,2)      -0.0218
   (3,2)       0.0873
   (4,2)      -0.0188
   (1,3)      -0.1528
   (2,3)      -0.0754
   (3,3)      -0.0899
   (4,3)       0.3181

sum(fit(A).^2)

ans = 0.0404

function res = fit(A)
  m0 = [0.5 0 0 0 ]';
  T = [299,1307.99];
  [T,m] = ode45(@(t,m) myfunction(t,m,A), T, m0);

  experiment=[0,0.02589,0.0540,0.01821];
  res = [m(end,1)-experiment(1);m(end,2)-experiment(2);...
         m(end,3)-experiment(3);m(end,4)-experiment(4)];
end

function dmdt = myfunction(T,m,A)
  R = 8.314*(10^(-3));
  E1 = 140;
  E2 = 133;
  E3 = 121;

  dmdt = zeros(4,1);
  dmdt(1) = -m(1)*(A(1)*exp(-E1/(R*T))+A(2)*exp(-E2/(R*T))+A(3)*exp(-E3/(R*T)));
  dmdt(2) = m(1)*A(1)*exp(-E1/(R*T));
  dmdt(3) = A(2)*exp(-E2/(R*T))* m(1);
  dmdt(4) = A(3)*exp(-E3/(R*T))* m(1);
end

[Rosalind]

I need to have the T written as T = 299 + (0.25*(1:1900)-0.25); to be able to change the timing and so on. when I put that instead of T = [299,1307.99]; I get the same A as A0.

Even when I use the code above and only change the initial mass i get the same A and A0.

[Torsten]

I need to have the T written as T = 299 + (0.25*(1:1900)-0.25); to be able to change the timing and so on.

Whether you write T = 299 + (0.25*(1:1900)-0.25) or T = [299,773.75] only changes the number of outputs from ODE45. Since you only need the last output for time = 773.75, it does not make a difference. Especially it doesn't make any changes in "timing and so on".

Note that the T in the list of inputs to myfunction is not the vector T = 299 + (0.25*(1:1900)-0.25), but only the actual integration time (thus a scalar between 299 and 773.75).

Even when I use the code above and only change the initial mass i get the same A and A0.

The MATLAB code works - it's up to you now to produce senseful results. I mentionned the problems with your data base.

[Rosalind]

the thing is that in order ro have the mass balance correct i should have put m0 = [0.098 0 0 0 ] instead of m0 = [0.5 0 0 0 ]. this was my mistake in the first place.And when i do change it I have the equal A0 and A, which is not correct.

[Torsten]

A0=1e-4*[1.5*(10^7),1*(10^6),1*(10^6)]

A0 = 1×3

1500 100 100

sum(fit(A0).^2)

ans = 0.0059

options = optimset('TolFun',1e-12,'TolX',1e-12);
[A,resnorm,residual,exitflag,output,lambda,Jacobian] = lsqnonlin(@fit,A0)

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

A = 1×3

1.0e+03 * 1.8842 1.1707 0.3470

resnorm = 0.0013

residual = 4×1

0.0225 -0.0044 -0.0270 0.0088

exitflag = 1

output = struct with fields:
    firstorderopt: 7.4758e-07
       iterations: 7
        funcCount: 32
     cgiterations: 0
        algorithm: 'trust-region-reflective'
         stepsize: 639.9974
          message: 'Local minimum found.↵↵Optimization completed because the size of the gradient is less than↵the value of the optimality tolerance.↵↵<stopping criteria details>↵↵Optimization completed: The first-order optimality measure, 7.475833e-07,↵is less than options.OptimalityTolerance = 1.000000e-06.'

lambda = struct with fields:
    lower: [3×1 double]
    upper: [3×1 double]

Jacobian = 
   1.0e-04 *

   (1,1)      -0.0511
   (2,1)       0.0952
   (3,1)      -0.0227
   (4,1)      -0.0214
   (1,2)      -0.1018
   (2,2)      -0.0382
   (3,2)       0.1841
   (4,2)      -0.0441
   (1,3)      -0.3331
   (2,3)      -0.1322
   (3,3)      -0.1615
   (4,3)       0.6267

sum(fit(A).^2)

ans = 0.0013

function res = fit(A)
  m0 = [0.098 0 0 0 ]';
  T = [299,1307.99];
  options = odeset('RelTol',1e-8,'AbsTol',1e-8);
  [T,m] = ode45(@(t,m) myfunction(t,m,A), T, m0, options);

  experiment=[0,0.02589,0.0540,0.01821];
  res = [m(end,1)-experiment(1);m(end,2)-experiment(2);...
         m(end,3)-experiment(3);m(end,4)-experiment(4)];
end

function dmdt = myfunction(T,m,A)
  R = 8.314*(10^(-3));
  E1 = 140;
  E2 = 133;
  E3 = 121;

  dmdt = zeros(4,1);
  dmdt(1) = -m(1)*(A(1)*exp(-E1/(R*T))+A(2)*exp(-E2/(R*T))+A(3)*exp(-E3/(R*T)));
  dmdt(2) = m(1)*A(1)*exp(-E1/(R*T));
  dmdt(3) = A(2)*exp(-E2/(R*T))* m(1);
  dmdt(4) = A(3)*exp(-E3/(R*T))* m(1);
end

[Rosalind]

Thank you Torsten. it worked.

Can you please tell me what to be careful about (or consider) if I want to expand my ODE equations and have more parameters to estimate?

[Torsten]

Can you please tell me what to be careful about (or consider) if I want to expand my ODE equations and have more parameters to estimate?

Try to get more experimental data (especially not only for the end time, but for the complete course of the process).