Fitting data with multiple inputs, ODE equation, and lsqnonlin

Hi everyone,

I am currently trying out MATLAB fitting to model a series chemical kinetics, which is based on an existing question in this forum seen here.

In my case, I add parameter of temperature to determine the reaction constant ki. The problem is intended to be solved using lsqnonlin.

Thus, the objective is to fit the data and parameters to find the values of Arrhenius coefficients A1-4 and Ea1-4.

Currently, I am having a code (seen at the end of this post) that generates some errors and can't be executed, so there I must've done some mistakes in either coding or my way of thinking. I looked up and refer to similar problems in the forum as well, but it seems it hasn't been working for my case so far.

https://www.mathworks.com/matlabcentral/answers/43439-monod-kinetics-and-curve-fitting

https://www.mathworks.com/matlabcentral/answers/592423-fitting-data-with-ode-equation-and-multiple-inputs?s_tid=answers_rc1-3_p3_Topic

I'm still pretty new in Matlab, especially in fitting nonlinear models. I hope someone can help me find the mistakes and give a solution to the problem.

Thanks in advance.

Warning: Length of lower bounds is < length(x); filling in missing lower bounds with -Inf. 
> In checkbounds (line 33)
In lsqnsetup (line 77)
In lsqnonlin (line 207)
In Untitled2_trial3 (line 53) 
Warning: Length of upper bounds is > length(x); ignoring extra bounds. 
> In checkbounds (line 41)
In lsqnsetup (line 77)
In lsqnonlin (line 207)
In Untitled2_trial3 (line 53) 
Exiting due to infeasibility: at least one lower bound exceeds the corresponding upper bound.
>> 

%Reaction system
%c1 --> c2 
%c1 --> p1
%c2 --> p2
%Reaction rates are described as r = k*c^n
%As c2 --> c3 is found negligable, data is available only for c1 and c2
%The data for p1 and p2 is unknown.

%Data
time = [10 20 40 60 100 140 180];
c1 = [0.508 0.442 0.385 0.354 0.299 0.267 0.258];
c2 = [0.246 0.301 0.359 0.398 0.422 0.467 0.489];
T_o = [400 423 433 456 467 468 463];
c1_0 = 0.711;
c2_0 = 0;
T_o_0 = 400;
R = 8.314;
%Plot
time_plot = [0,time];
c1_plot = [c1_0,c1];
c2_plot = [c2_0,c2];
T_o_plot = [0,T_o];

%figure 1, C1 and C2 vs time
figure(1)
plot(time_plot,c1_plot,'bo')
grid on
hold on
plot(time_plot,c2_plot,'ro')
hold off
xlabel('Time [min]')
ylabel('Concentration [mol/L]')
legend('c1','c2')

%figure 2, C1 and C2 vs T
figure(2)
plot(T_o_plot,c1_plot,'bo')
grid on
hold on
plot(T_o_plot,c2_plot,'ro')
hold off
xlabel('Temperature [K]')
ylabel('Concentration [mol/L]')
legend('c1','c2')

%Data processing
x_in = [time',T_o'];
y_in = [c1',c2'];

%Lsqnonlin
%[p,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqnonlin(@diff1,p0,x_in,y_in,lb,ub);
%FitData = diff1(p,time,1);
%c1_out = FitData(:,1);
%c2_out = FitData(:,2);

% Options for fitting
options = optimoptions('lsqnonlin', ...
    'Display', 'iter', ...
    'Algorithm', 'levenberg-marquardt', ...
    'UseParallel', false, ...
    'StepTolerance', 1e-6);
%initial value for iteration
B0 = [0.1,1,0.1,1,0.1,1,0.1,1,0.1,1,0.1,1,c1(1),c2(1),T_o(1)];


% lsqnonlin
optimizedParams = lsqnonlin(@diff1, B0, x_in, y_in, [], [], options);

function C = diff1(p,time,~)
%dc1/dt = -(k1*c1^n1+k3*c1^k3)
%dc2/dt = k1*c1^n1-(k2*c2^k2+k4*c2^n4)
%where
%k1 = A1*exp(-E1/RT_o)
%k2 = A2*exp(-E2/RT_o)
%k3 = A3*exp(-E3/RT_o)
%k4 = A4*exp(-E4/RT_o)
%
%variables: c1 = x(1), c2 = x(2), T_o = x(3)
%parameters: 
%A1 = B(1), E1 = B(2), n1 = B(3)
%A2 = B(4), E2 = B(5), n2 = B(6)
%A3 = B(7), E3 = B(8), n3 = B(9)
%A4 = B(10), E4 = B(11), n4 = B(12)


x0 = B(11:12);
[t,Cout] = ode45(@DifEq, time, x0);
      function dC = DifEq(t, x)
          xdot = zeros(3,1);
          xdot(1) = -(B(1)*exp(-B(2)/R/x(3)).*x(1).^B(3)-B(7)*exp(-B(8)/R/x(3)).*x(1)^B(9));
          xdot(2) = B(4)*exp(-B(5)/R/x(3)).*x(1).^B(6)-B(10)*exp(-B(11)/R/x(3)).*x(2).^B(12);        
          dC = xdot;
      end
C = Cout;
end

Answer 1

[Star Strider]

The fundamental problem is that lsqnonlin ~= lsqcurvefit! They are definitely not interchangable, and have different argument lists and calling conventions, the ones used here being for lsqcurvefit, so that change was straightforward (after checking to be sure everything matched correctly).

Beyond that, there is no ‘x(3)’ since there is no third differential equation. Since ‘x(3)’ should be temperature (that I here assume is a funciton of time), I created:

Temp = interp1(x_in(:,1), x_in(:,2),t);                                               % <— ADDED & Changed ‘x(3)’ To ‘Temp’

to interpolate temperature with time, and then made the substitutions in the differential equations. If temperature is not a function of time, it will be necessary to create a function that expresses temperature in a way that is meaningful for the differential equations.

I also set the lower bounds on the parameters to zero, since concentrations and kinetic constants, at least in this system, appear to always be greater than zero.

This version runs without errors:

%Reaction system
%c1 --> c2 
%c1 --> p1
%c2 --> p2
%Reaction rates are described as r = k*c^n
%As c2 --> c3 is found negligable, data is available only for c1 and c2
%The data for p1 and p2 is unknown.

%Data
time = [10 20 40 60 100 140 180];
c1 = [0.508 0.442 0.385 0.354 0.299 0.267 0.258];
c2 = [0.246 0.301 0.359 0.398 0.422 0.467 0.489];
T_o = [400 423 433 456 467 468 463];
c1_0 = 0.711;
c2_0 = 0;
T_o_0 = 400;
R = 8.314;

%Data processing
x_in = [time(:),T_o(:)];
y_in = [c1(:),c2(:)];

%Lsqnonlin
%[p,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqnonlin(@diff1,p0,x_in,y_in,lb,ub);
%FitData = diff1(p,time,1);
%c1_out = FitData(:,1);
%c2_out = FitData(:,2);

% Options for fitting
options = optimoptions('lsqcurvefit', ...
    'Display', 'iter', ...
    'Algorithm', 'levenberg-marquardt', ...
    'UseParallel', false, ...
    'StepTolerance', 1e-6);
%initial value for iteration
B0 = [0.1,1,0.1,1,0.1,1,0.1,1,0.1,1,0.1,1,c1(1),c2(1),T_o(1)];


% lsqnonlin
optimizedParams = lsqcurvefit(@diff1, B0, x_in, y_in, zeros(size(B0)), [], options)

function C = diff1(B,x_in)
%dc1/dt = -(k1*c1^n1+k3*c1^k3)
%dc2/dt = k1*c1^n1-(k2*c2^k2+k4*c2^n4)
%where
%k1 = A1*exp(-E1/RT_o)
%k2 = A2*exp(-E2/RT_o)
%k3 = A3*exp(-E3/RT_o)
%k4 = A4*exp(-E4/RT_o)
%
%variables: c1 = x(1), c2 = x(2), T_o = x(3)
%parameters: 
%A1 = B(1), E1 = B(2), n1 = B(3)
%A2 = B(4), E2 = B(5), n2 = B(6)
%A3 = B(7), E3 = B(8), n3 = B(9)
%A4 = B(10), E4 = B(11), n4 = B(12)

time = x_in(:,1);

x0 = B(11:12);
[t,Cout] = ode45(@DifEq, time, x0);
      function dC = DifEq(t, x)
          xdot = zeros(2,1);
          Temp = interp1(x_in(:,1), x_in(:,2),t);                                               % <— ADDED & Changed 'x(3)' To 'Temp'
          xdot(1) = -(B(1)*exp(-B(2)/R/Temp).*x(1).^B(3)-B(7)*exp(-B(8)/R/Temp).*x(1)^B(9));
          xdot(2) = B(4)*exp(-B(5)/R/Temp).*x(1).^B(6)-B(10)*exp(-B(11)/R/Temp).*x(2).^B(12);
          dC = xdot;
      end
C = Cout;
end

%Plot
time_plot = [time];
c1_plot = [c1];
c2_plot = [c2];
T_o_plot = [T_o];

Cest = diff1(optimizedParams,x_in)

%figure 1, C1 and C2 vs time
figure(1)
plot(time_plot,c1_plot,'bo')
grid on
hold on
plot(time_plot,c2_plot,'ro')
plot(x_in(:,1), Cest, '--')
hold off
xlabel('Time [min]')
ylabel('Concentration [mol/L]')
legend('c1','c2')

%figure 2, C1 and C2 vs T
figure(2)
plot(T_o_plot,c1_plot,'bo')
grid on
hold on
plot(T_o_plot,c2_plot,'ro')
plot(x_in(:,2), Cest, '--')
hold off
xlabel('Temperature [K]')
ylabel('Concentration [mol/L]')
legend('c1','c2')

I shifted the plots to the end in order to incorporate the fitted equations, and added those.

It will be necessary to experiment with different initial parameter estimates to get a decent fit, since the current estimates do not create that. If you have the Golbal Optimization Toolbox, use the ga (genetic algorithm) function for this. It will generally find a reasonable fit, although the 'InitialPopulationMatrix' may need to be designed specifically for these parameters. I can hel with that if necessary.

[Anthony Hamzah]

Thanks for the helpful response, much appreciated. I'm currently trying out different initial estimates and learning GA toolbox in the same time. I will reach you out later.

Just a simple favor, I have difficulty displaying the fitted constant values of B (B(1) to (12)). I tried by typing

disp(B)

but somehow it didn't work at all (the variable is deemed unrecognized). Could you tell me how to display the value?

[Star Strider]

My pleasure!

If you want the estimated parameters:

fprintf('Parameters: \n')
for k = 1:numel(optimizedParams)
    fprintf('\tB(%2d) = %23.15E\n', k, optimizedParams(k))
end

will print them to theCommand Window.

This produces:

Parameters: 
	B( 1) =   9.997794030302047E-02
	B( 2) =   1.000000075249041E+00
	B( 3) =   1.007603515365702E-01
	B( 4) =   1.000025114954299E+00
	B( 5) =   9.999999564401417E-02
	B( 6) =   9.999757769443918E-01
	B( 7) =   1.024040170625801E-01
	B( 8) =   9.999999935095988E-01
	B( 9) =   9.978242903083580E-02
	B(10) =   9.999846654525757E-01
	B(11) =   1.003304243186730E-01
	B(12) =   9.999484385750329E-01
	B(13) =   5.080000000000000E-01
	B(14) =   2.460000000000000E-01
	B(15) =   4.000000000000000E+02

when I just now ran it.

I have prototype ga code for these sorts of problems, so if you’re interested, I’ll adapt it to this problem and post it. It only makes sense to do that if you have the Global Optimization Toolbox, so that you can run it and experiment with it.

The exponents should likely be integers (or so I understood from my undergraduate physical chemistry courses), and that can be forced by rounding them inside the objective function. I need to know if you want to do that, and what ‘B’ elements are the exponents (so I don’t have to dissect your code to figure out what they are).

[Anthony Hamzah]

Thanks a lot.

I've been experimenting with various initial guess, but it seems none fits the data so far. The iterated values constantly give numbers that are way too close to the initial guess.

On the other hand, it turns out the Global Optimization toolbox is installed in my license. Would be very happy if you post the example of the GA code you mentioned.

[Star Strider]

My pleasure!

This is the code I’m curently running:

function Anthony_Hamzah_GA
% 2016 12 03
% NOTES: 
% 
%   1. The ‘theta’ (parameter) argument has to be first in your 
%       ‘kinetics’ function,
%   2. You need to return ALL the values from ‘DifEq’ since you are fitting
%       all the values

function C = diff1(B,x_in)
%dc1/dt = -(k1*c1^n1+k3*c1^k3)
%dc2/dt = k1*c1^n1-(k2*c2^k2+k4*c2^n4)
%where
%k1 = A1*exp(-E1/RT_o)
%k2 = A2*exp(-E2/RT_o)
%k3 = A3*exp(-E3/RT_o)
%k4 = A4*exp(-E4/RT_o)
%
%variables: c1 = x(1), c2 = x(2), T_o = x(3)
%parameters: 
%A1 = B(1), E1 = B(2), n1 = B(3)
%A2 = B(4), E2 = B(5), n2 = B(6)
%A3 = B(7), E3 = B(8), n3 = B(9)
%A4 = B(10), E4 = B(11), n4 = B(12)
time = x_in(:,1);
x0 = B(11:12);
[t,Cout] = ode45(@DifEq, time, x0);
      function dC = DifEq(t, x)
          xdot = zeros(2,1);
          Temp = interp1(x_in(:,1), x_in(:,2),t);                                               % <— ADDED & Changed 'x(3)' To 'Temp'
          xdot(1) = -(B(1)*exp(-B(2)/R/Temp).*x(1).^B(3)-B(7)*exp(-B(8)/R/Temp).*x(1)^B(9));
          xdot(2) = B(4)*exp(-B(5)/R/Temp).*x(1).^B(6)-B(10)*exp(-B(11)/R/Temp).*x(2).^B(12);
          dC = xdot;
      end
C = Cout;
end



%Reaction system
%c1 --> c2 
%c1 --> p1
%c2 --> p2
%Reaction rates are described as r = k*c^n
%As c2 --> c3 is found negligable, data is available only for c1 and c2
%The data for p1 and p2 is unknown.
%Data
time = [10 20 40 60 100 140 180];
c1 = [0.508 0.442 0.385 0.354 0.299 0.267 0.258];
c2 = [0.246 0.301 0.359 0.398 0.422 0.467 0.489];
T_o = [400 423 433 456 467 468 463];
c1_0 = 0.711;
c2_0 = 0;
T_o_0 = 400;
R = 8.314;
%Data processing
x_in = [time(:),T_o(:)];
y_in = [c1(:),c2(:)];

% [theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
ftns = @(B) norm(y_in-diff1(B,x_in));
PopSz = 500;
Parms = 15;
opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',2E3, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);
t0 = clock;
fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)
[B,fval,exitflag,output] = ga(ftns, Parms, [],[],[],[],zeros(Parms,1),Inf(Parms,1),[],[],opts)
t1 = clock;
fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)
GA_Time = etime(t1,t0)
DT_GA_Time = datetime([0 0 0 0 0 GA_Time], 'Format','HH:mm:ss.SSS');
fprintf('\nElapsed Time: %23.15E\t\t%s\n\n', GA_Time, DT_GA_Time)

fprintf(1,'\tParameters:\n')
for k1 = 1:numel(B)
    fprintf(1, '\t\tB(%2d) = %11.8f\n', k1, B(k1))
end


%Plot
time_plot = [time];
c1_plot = [c1];
c2_plot = [c2];
T_o_plot = [T_o];

Cest = diff1(optimizedParams,x_in)

%figure 1, C1 and C2 vs time
figure(1)
plot(time_plot,c1_plot,'bo')
grid on
hold on
plot(time_plot,c2_plot,'ro')
plot(x_in(:,1), Cest, '--')
hold off
xlabel('Time [min]')
ylabel('Concentration [mol/L]')
legend('c1','c2')

%figure 2, C1 and C2 vs T
figure(2)
plot(T_o_plot,c1_plot,'bo')
grid on
hold on
plot(T_o_plot,c2_plot,'ro')
plot(x_in(:,2), Cest, '--')
hold off
xlabel('Temperature [K]')
ylabel('Concentration [mol/L]')
legend('c1','c2')

end

It appears to work correctly, however it is still in the first generation so I have no idea what it will do or how long it will take to converge. Optimising 15 parameters is not trivial.

It’s the end of my day here, so I’ll let run all night and see what it does. Windows 10 on this computer has the unfortunate habit of crashing unpredictably, so unless that occurs, I should have results tomorrow morning. I definitely encourage you to run it as well unless you need to use MATLAB for something else for the time being.

(The edit corrected some typographical errors that I didn’t initially see.)

[Anthony Hamzah]

Once again, thanks a lot.

Tried to run the GA simulation, it turns out the 'opts =...' is somehow erroneous. Still don't know why it doesn't work.

Error in Anthony_Hamzah_GA (line 55)
opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3,
'MaxGenerations',2E3, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);

[Star Strider]

It is taking forever to converge, however it is running without error.

The optimoptions call is not throwing any errors when I run my code.

I found an error in the original code that I didn’t see previously, so ‘x0’ should be:

x0 = B(13:14);

and there are only 14 parameters as the result.

Starting it over with these changes.

... Later that same day ...

EDIT — (25 Nov 2020 at 16:54)

I made a few more changes and finally got it running and giving reasonable results. Note that the ‘PopSize’ is reduced and the number of generations increased. This is necessary because of the number of parameters. There are also changes in the optimoptions call to provide reaonably appropriate initial magnitudes for the parameters, and to decrease the stall time. The ‘x_in_ext’ matrix increases the precision of the fitted curve in the plots.

The current code (including an additional figure(3) presenting both original plots in one figure window):

function Anthony_Hamzah_GA
% 2016 12 03
% NOTES: 
% 
%   1. The ‘theta’ (parameter) argument has to be first in your 
%       ‘kinetics’ function,
%   2. You need to return ALL the values from ‘DifEq’ since you are fitting
%       all the values

function C = diff1(B,x_in)
%dc1/dt = -(k1*c1^n1+k3*c1^k3)
%dc2/dt = k1*c1^n1-(k2*c2^k2+k4*c2^n4)
%where
%k1 = A1*exp(-E1/RT_o)
%k2 = A2*exp(-E2/RT_o)
%k3 = A3*exp(-E3/RT_o)
%k4 = A4*exp(-E4/RT_o)
%
%variables: c1 = x(1), c2 = x(2), T_o = x(3)
%parameters: 
%A1 = B(1), E1 = B(2), n1 = B(3)
%A2 = B(4), E2 = B(5), n2 = B(6)
%A3 = B(7), E3 = B(8), n3 = B(9)
%A4 = B(10), E4 = B(11), n4 = B(12)
time = x_in(:,1);
x0 = B(13:14);
[t,Cout] = ode45(@DifEq, time, x0);
      function dC = DifEq(t, x)
          xdot = zeros(2,1);
          Temp = interp1(x_in(:,1), x_in(:,2),t);                                               % <— ADDED & Changed 'x(3)' To 'Temp'
          xdot(1) = -(B(1)*exp(-B(2)/R/Temp).*x(1).^B(3)-B(7)*exp(-B(8)/R/Temp).*x(1)^B(9));
          xdot(2) = B(4)*exp(-B(5)/R/Temp).*x(1).^B(6)-B(10)*exp(-B(11)/R/Temp).*x(2).^B(12);
          dC = xdot;
      end
C = Cout;
end



%Reaction system
%c1 --> c2 
%c1 --> p1
%c2 --> p2
%Reaction rates are described as r = k*c^n
%As c2 --> c3 is found negligable, data is available only for c1 and c2
%The data for p1 and p2 is unknown.
%Data
time = [10 20 40 60 100 140 180];
c1 = [0.508 0.442 0.385 0.354 0.299 0.267 0.258];
c2 = [0.246 0.301 0.359 0.398 0.422 0.467 0.489];
T_o = [400 423 433 456 467 468 463];
c1_0 = 0.711;
c2_0 = 0;
T_o_0 = 400;
R = 8.314;
%Data processing
x_in = [time(:),T_o(:)];
y_in = [c1(:),c2(:)];

% [theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
ftns = @(B) norm(y_in-diff1(B,x_in));
PopSz = 10;
Parms = 14;
opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms).*[1E-05, 1E-04, 1E-05, 1E-04, 1E-05, 1E-04, 1E-05, 1E-04, 1E-05, 1E-04, 1E-05, 1E-04, 1E-04, 1E-05]*0.1, 'MaxGenerations',5E3, 'PlotFcn',@gaplotbestf, 'PlotInterval',1, 'StallTimeLimit',60);
t0 = clock;
fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)
[B,fval,exitflag,output] = ga(ftns, Parms, [],[],[],[],zeros(Parms,1),Inf(Parms,1),[],[],opts)
t1 = clock;
fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)
GA_Time = etime(t1,t0)
DT_GA_Time = datetime([0 0 0 0 0 GA_Time], 'Format','HH:mm:ss.SSS');
fprintf('\nElapsed Time: %23.15E\t\t%s\n\n', GA_Time, DT_GA_Time)

fprintf(1,'\tParameters:\n')
for k1 = 1:length(B)
    fprintf(1, '\t\tB(%2d) = %11.8f\n', k1, B(k1))
end

%Plot
time_plot = [time];
c1_plot = [c1];
c2_plot = [c2];
T_o_plot = [T_o];

N = 50;
timev = linspace(min(time), max(time), N);
Tempv = linspace(min(T_o), max(T_o), N);
x_in_ext = [timev(:), Tempv(:)];
Cest = diff1(B,x_in_ext);

%figure 1, C1 and C2 vs time
figure(1)
plot(time_plot,c1_plot,'bo')
grid on
hold on
plot(time_plot,c2_plot,'ro')
plot(x_in_ext(:,1), Cest(:,1), '-b')
plot(x_in_ext(:,1), Cest(:,2), '-r')
hold off
xlabel('Time [min]')
ylabel('Concentration [mol/L]')
legend('c1','c2')

%figure 2, C1 and C2 vs T
figure(2)
plot(T_o_plot,c1_plot,'bo')
grid on
hold on
plot(T_o_plot,c2_plot,'ro')
plot(x_in_ext(:,2), Cest(:,1), '-b')
plot(x_in_ext(:,2), Cest(:,2), '-r')
hold off
xlabel('Temperature [K]')
ylabel('Concentration [mol/L]')
legend('c1','c2')

figure(3)
subplot(2,1,1)
plot(time_plot,c1_plot,'bo')
grid on
hold on
plot(time_plot,c2_plot,'ro')
plot(x_in_ext(:,1), Cest(:,1), '-b')
plot(x_in_ext(:,1), Cest(:,2), '-r')
hold off
xlabel('Time [min]')
ylabel('Concentration [mol/L]')
legend('c1','c2')

%figure 2, C1 and C2 vs T
subplot(2,1,2)
plot(T_o_plot,c1_plot,'bo')
grid on
hold on
plot(T_o_plot,c2_plot,'ro')
plot(x_in_ext(:,2), Cest(:,1), '-b')
plot(x_in_ext(:,2), Cest(:,2), '-r')
hold off
xlabel('Temperature [K]')
ylabel('Concentration [mol/L]')
legend('c1','c2')

end

producing thesse parameters:

	Parameters:
		B( 1) =  0.51277329
		B( 2) =  0.08625777
		B( 3) =  0.01908611
		B( 4) =  0.03112946
		B( 5) =  0.20644214
		B( 6) =  0.59838085
		B( 7) =  0.51375554
		B( 8) =  0.40580142
		B( 9) =  0.02256135
		B(10) =  0.01967728
		B(11) =  0.06631409
		B(12) =  0.31020972
		B(13) =  0.40685799
		B(14) =  0.28383620

with this fitness value (norm of the residuals):

fval =
   120.8957e-003

and this plot for figure(3):

This is typical for the best results I’m getting, so likely the best the algorithm can do. The only suggestion I have is to review the differential equation system to be certain they are coded correctly. I did not force the order exponents to be integers, since I am not certain which ones they are. It might be worthwhile to include ‘c3’ since it could further define ‘k2’, ‘c2’ and ‘n2’.

.

[Anthony Hamzah]

Glad it works in your computer, at least I know what to expect once the problem is sorted out. However it still generates similar error with the 'optimoptions' syntax to that I reported previously

Anthony_Hamzah_GA1
Error using optimoptions (line 124)
Invalid solver specified. Provide a solver name or handle (such as 'fmincon' or @fminunc).
Type DOC OPTIMOPTIONS for a list of solvers.

Error in Anthony_Hamzah_GA1 (line 54)
opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms).*[1E-05, 1E-04,
1E-05, 1E-04, 1E-05, 1E-04, 1E-05, 1E-04, 1E-05, 1E-04, 1E-05, 1E-04, 1E-04, 1E-05]*0.1, 'MaxGenerations',5E3,
'PlotFcn',@gaplotbestf, 'PlotInterval',1, 'StallTimeLimit',60);

I honestly have no idea why it works in one computer and not in the other one despite both having Global Optimization Toolbox installed.

EDIT — (26 Nov 2020 at JST 10:36)

Finally able to make it work by uninstalling the tollbox and later reinstalling it. The result is rather different though despite the same population quantity and I'm increasing the population size to make it closer to the concentration data.

where B values

Stop Time: 2020-11-26 10:36:08.8080

GA_Time =

   60.7290


Elapsed Time:   6.072900000000000E+01		00:01:00.729

	Parameters:
		B( 1) =  0.01420597
		B( 2) =  0.14147360
		B( 3) =  0.02138446
		B( 4) =  0.09937078
		B( 5) =  0.19028512
		B( 6) =  0.05423926
		B( 7) =  0.01956314
		B( 8) =  0.13490452
		B( 9) =  0.45489825
		B(10) =  0.11626381
		B(11) =  0.04371411
		B(12) =  0.24883097
		B(13) =  0.46448139
		B(14) =  0.42516848

[Star Strider]

The ga function will not always converge to the same results, especially with more difficult problems such as this. It will be necessary to run it perhaps several times, and choose the result with the lowest fitness, that (in my code) will produce the best fit to the data. Note that my best result produced a fitness value of about 0.121.

[Anthony Hamzah]

After running the code with population size of 20 several times, this is the best I can get.

By the way, could you show me how to print the fitness and R-square values?

[Star Strider]

The fitness value is ‘fval’ from the ga function output.

The value would be:

RMSE = sqrt(mean(y_in-diff1(B,x_in).^2));
OLSCF = @(B) sum(y_in-diff1(B,x_in)).^2);
SStot = sum((y - mean(y)).^2);
SSres = OLSCF(B);
Rsq   = 1 - (SSres/SStot);                      % Compute R-squared

I did not specifically compute this with my code. It should work.

[Anthony Hamzah]

Sorry for the late reply. I've been busy handling other stuff at work all day yesterday.

Thanks for the help, I find it extremely helpful.

[Star Strider]

No worries!

As always, my pleasure!

Answer 2

[Alan Weiss]

I do not have the time right now to help debug your code, sorry. But I believe that you might be able to use some documentation examples that run correctly as models to help figure out your own error:

Fit ODE, Problem-Based

Fit an Ordinary Differential Equation (ODE)

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation