How to fit kinetic model to estimate kinetic parameter from experimental data

I am getting errors in the output. Please help be solve the problem. I have attached the Excel file of the experimental data.

Thank you.

%Fermentation data

Xdata = readmatrix('Batch1.xlsx','Sheet','sheet1','Range','e2:e13');
Gdata = readmatrix('Batch1.xlsx','Sheet','sheet1','Range','b2:b13');
Bdata = readmatrix('Batch1.xlsx','Sheet','sheet1','Range','c2:c13');
Edata = readmatrix('Batch1.xlsx','Sheet','sheet1','Range','d2:d13');


% Group all data into yariable y
yinv =[Xdata'; Gdata'; Bdata'; Edata'];


%Data for time
timeex = readmatrix('batch1.xlsx','Sheet','sheet1','Range','a2:a13');

%Set ODE solver
% Initial cell biomass concentration (Initial condition)
y0=[1.04 0.78 0.00 0.00];

%Fermentation time
tspan = linspace(0,72);

%Set optimization problem
%Let b = matrix of paramters to be determined
% b= [um ks k1 K1G K1B k2 K2G YXG m k3]
b = optimvar('b',10,"LowerBound",0,"UpperBound",72);

%Set functions for ODE solver for solving ODE
function solferment = toODE(b,tspan,y0)
    sol = ode45(@(t,y)batchferment(t,y,b),tspan,y0);
    solferment = deval(sol,tspan);
end

%Convert function for ODE solving to optimization expression 

%To use RtoODE in an objective function, convert the function to an
%optimization expression by using fcn2optimexpr. 
myfcn = fcn2optimexpr(@toODE,b,timeex,y0);

%Express the objective function as the sum of squared differences between
%the ODE solution and the solution with true parameters.
SSE = sum(sum((myfcn-yinv).^2));

%Create an optimization problem with the objective function SSE.
prob = optimproblem ("Description","Fit ODE parameters",'ObjectiveSense','min');

%Objective function (to be minimized)
prob.Objective = SSE;

%Show structure of problem
show(prob)

%Solve Problem 
%To find the best-fitting parameters x, give an initial guess
%x0 for the solver and call solve. 
% Set initial guesses of parameters
initialGuess.b = [0.18 1.0 0.61 0.18 5.85 3.20 16.25 0.11 3.40 3.02];

%Solve optimization problem
[sol,optval] = solve(prob,initialGuess);

%Extract the results
%Fitted Parameters
bfinal =sol.b;

%Sum of square Error
SSEfinal = optval;


%Plot the simulated data and experimental data
%Call ODE to solve an equation using Final Fitted Parameters (bfinal)
solysim = ode45(@(t,y)batchferment(t,y,bfinal),tspan,y0);
%Evaluate the simulated results at specified tspan
ysim = deval(solysim,tspan);

%Plot graphs
figure;
plot(timeex,Xdata,'*r',tspan,ysim(1,:),'-b')
legend('exp','sim')
xlabel('Time (h)')
ylabel('X (g/L)')

figure;
plot(timeex,Gdata','*r',tspan,ysim(2,:),'-b')
legend('exp','sim')
xlabel('Time (h)')
ylabel('G (g/L)')

figure;
plot(timeex,Bdata','*r',tspan,ysim(3,:),'-b')
legend('exp','sim')
xlabel('Time (h)')
ylabel('B (g/L)')

figure;
plot(timeex,Edata','*r',tspan,ysim(4,:),'-b')
legend('exp','sim')
xlabel('Time (h)')
ylabel('E (g/L)')


%Equations for Batch 
%y(1) = X = Biomass Concentration (g/l)
%y(2) = G = Glucose Concentration (g/l)
%y(3) = B = Cellobiose Concentration (g/L)
%y(4) = E = Ethanol Concentration (g/l)

function dydt = batchferment(t,y,b)

   y(1) = X; 
   y(2) = G;
   y(3) = B;   
   y(4) = E;

    %%%Growth equations
    %dx/dt = (um*X*G)/(KG+G)
    u = (b(1)*y(1)*y(2))/(b(2)+y(2));
    %u = (b(1)*X*G)/(b(2)+G);
    %parameters
  b(1) = um (1/h);
  b(2) = ks (g/L);


    %%%Cellobiose equations
    C = b(11)-(0.9*y(2))-(0.947*y(3))-(0.9*y(4)/0.511)-(1.137*(y(1)-1.04));
     r1 = (b(3)*C)/(1+(y(2)/b(4))+(y(3)/b(5)));
     r2 = (b(6)*y(3))/(1+(y(2)/b(7)));
  b(3) = k1;
  b(4) = K1G;
  b(5) = K1B;
  b(6) = k2;
  b(7) = K2G;

  %C = Cellulose concentration (g/l)
  %C = C0-(0.9*G)-(0.947*B)-(0.9*E/0.511)-(1.137*(X-X0))
  %r1 = (k1*C)/(1+(G/K1G)+(B/K1B))
  %r2 = (k2*B)/(1+(G/K2G)

 %C = C0-(0.9*G)-(0.947*B)-(0.9*E/0.511)-(1.137*(X-X0))
  %C = 83.65-(0.9*y(2))-(0.947*y(3))-(0.9*y(4)/0.511)-(1.137*(y(1)-1.04))
  
    
    %%% Glucose equations
    rG = ((1/b(8))*(dx/dt))-(b(9)*y(1));
  b(8) = YXG;
  b(9) = m;
  %rG = ((1/YXG)*(dX/dt))-(m*X)

    
    %%% Ethanol Concentration
    Et = k3*(u*X)*(1/YXG);
  b(10) = k3;
  

    %Material balance equations
    dXdt = u*X;
    dGdt = (r2/0.95)-rG;
    dBdt = (r1/0.947)-r2;
    dEdt = Et;
    dydt = [dXdt;dGdt;dBdt;dEdt];
end

This is the output after running:

>> Batch1

Unrecognized function or variable 'X'.

Error in Batch1>batchferment (line 110)

y(1) = X;

Error in Batch1>@(t,y)batchferment(t,y,b) (line 30)

sol = ode45(@(t,y)batchferment(t,y,b),tspan,y0);

Error in odearguments (line 92)

f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 104)

odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error in Batch1>toODE (line 30)

sol = ode45(@(t,y)batchferment(t,y,b),tspan,y0);

Error in optim.problemdef.fcn2optimexpr

Error in optim.problemdef.fcn2optimexpr

Error in fcn2optimexpr

Error in Batch1 (line 38)

myfcn = fcn2optimexpr(@toODE,b,timeex,y0);

Caused by:

Function evaluation failed while attempting to determine output size. The function might contain an error, or might not be well-defined at

the automatically-chosen point. To specify output size without function evaluation, use 'OutputSize'.

>>

Answer 1

[Star Strider]

I got this to run (there were a number of coding errors, most of which I corrected). However it takes too long to run here, so there could be some coding errors I have not yet caught, although the rest of the code looks to me to be correct. I will try running it on MATLAB Online and post any further corrections (ir any are needed) as EDIT notes.

Try this —

%Fermentation data

Xdata = readmatrix('Batch1.xlsx','Sheet','sheet1','Range','e2:e13');
Gdata = readmatrix('Batch1.xlsx','Sheet','sheet1','Range','b2:b13');
Bdata = readmatrix('Batch1.xlsx','Sheet','sheet1','Range','c2:c13');
Edata = readmatrix('Batch1.xlsx','Sheet','sheet1','Range','d2:d13');


% Group all data into yariable y
yinv =[Xdata'; Gdata'; Bdata'; Edata'];


%Data for time
timeex = readmatrix('Batch1.xlsx','Sheet','sheet1','Range','a2:a13');

%Set ODE solver
% Initial cell biomass concentration (Initial condition)
y0=[1.04 0.78 0.00 0.00];

%Fermentation time
tspan = linspace(0,72);

%Set optimization problem
%Let b = matrix of paramters to be determined
% b= [um ks k1 K1G K1B k2 K2G YXG m k3]
b = optimvar('b',10,"LowerBound",0,"UpperBound",72);

%Set functions for ODE solver for solving ODE
function solferment = toODE(b,tspan,y0)
    sol = ode45(@(t,y)batchferment(t,y,b),tspan,y0);
    solferment = deval(sol,tspan);
end

%Convert function for ODE solving to optimization expression 

%To use RtoODE in an objective function, convert the function to an
%optimization expression by using fcn2optimexpr. 
myfcn = fcn2optimexpr(@toODE,b,timeex,y0);

%Express the objective function as the sum of squared differences between
%the ODE solution and the solution with true parameters.
SSE = sum(sum((myfcn-yinv).^2));

%Create an optimization problem with the objective function SSE.
prob = optimproblem ("Description","Fit ODE parameters",'ObjectiveSense','min');

%Objective function (to be minimized)
prob.Objective = SSE;

%Show structure of problem
show(prob)

%Solve Problem 
%To find the best-fitting parameters x, give an initial guess
%x0 for the solver and call solve. 
% Set initial guesses of parameters
initialGuess.b = [0.18 1.0 0.61 0.18 5.85 3.20 16.25 0.11 3.40 3.02];

%Solve optimization problem
[sol,optval] = solve(prob,initialGuess);

%Extract the results
%Fitted Parameters
bfinal =sol.b;

%Sum of square Error
SSEfinal = optval;


%Plot the simulated data and experimental data
%Call ODE to solve an equation using Final Fitted Parameters (bfinal)
solysim = ode45(@(t,y)batchferment(t,y,bfinal),tspan,y0);
%Evaluate the simulated results at specified tspan
ysim = deval(solysim,tspan);

%Plot graphs
figure;
plot(timeex,Xdata,'*r',tspan,ysim(1,:),'-b')
legend('exp','sim')
xlabel('Time (h)')
ylabel('X (g/L)')

figure;
plot(timeex,Gdata','*r',tspan,ysim(2,:),'-b')
legend('exp','sim')
xlabel('Time (h)')
ylabel('G (g/L)')

figure;
plot(timeex,Bdata','*r',tspan,ysim(3,:),'-b')
legend('exp','sim')
xlabel('Time (h)')
ylabel('B (g/L)')

figure;
plot(timeex,Edata','*r',tspan,ysim(4,:),'-b')
legend('exp','sim')
xlabel('Time (h)')
ylabel('E (g/L)')


%Equations for Batch 
%y(1) = X = Biomass Concentration (g/l)
%y(2) = G = Glucose Concentration (g/l)
%y(3) = B = Cellobiose Concentration (g/L)
%y(4) = E = Ethanol Concentration (g/l)

function dydt = batchferment(t,y,b)

   % y(1) = X; 
   % y(2) = G;
   % y(3) = B;   
   % y(4) = E;
   X = y(1);
   G = y(2);
   B = y(3);
   E = y(4);

    %%%Growth equations
    %dx/dt = (um*X*G)/(KG+G)
    u = (b(1)*y(1)*y(2))/(b(2)+y(2));
    %u = (b(1)*X*G)/(b(2)+G);
    %parameters
  % b(1) = um (1/h);
  % b(2) = ks (g/L);
  um = b(1);
  ks = b(2);


    %%%Cellobiose equations
    C = b(1)-(0.9*y(2))-(0.947*y(3))-(0.9*y(4)/0.511)-(1.137*(y(1)-1.04));
     r1 = (b(3)*C)/(1+(y(2)/b(4))+(y(3)/b(5)));
     r2 = (b(6)*y(3))/(1+(y(2)/b(7)));
  % b(3) = k1;
  % b(4) = K1G;
  % b(5) = K1B;
  % b(6) = k2;
  % b(7) = K2G;
  k1 = b(3);
  K1G = b(4);
  K1B = b(5);
  k2 = b(6);
  k2G = b(7);
  
  %C = Cellulose concentration (g/l)
  %C = C0-(0.9*G)-(0.947*B)-(0.9*E/0.511)-(1.137*(X-X0))
  %r1 = (k1*C)/(1+(G/K1G)+(B/K1B))
  %r2 = (k2*B)/(1+(G/K2G)

 %C = C0-(0.9*G)-(0.947*B)-(0.9*E/0.511)-(1.137*(X-X0))
  %C = 83.65-(0.9*y(2))-(0.947*y(3))-(0.9*y(4)/0.511)-(1.137*(y(1)-1.04))
  
    
    %%% Glucose equations
    % rG = ((1/b(8))*(dx/dt))-(b(9)*y(1));
  % b(8) = YXG;
  % b(9) = m;
  YXG = b(8)
  m = b(9);
  %rG = ((1/YXG)*(dX/dt))-(m*X)

    
    %%% Ethanol Concentration
    
    % b(10) = k3;
    k3 = b(10); 
    Et = k3*(u*X)*(1/YXG);  

    %Material balance equations
    dXdt = u*X;
    rG = ((1/b(8))*(dXdt))-(b(9)*y(1));
    dGdt = (r2/0.95)-rG;
    dBdt = (r1/0.947)-r2;
    dEdt = Et;
    dydt = [dXdt;dGdt;dBdt;dEdt];
end

EDIT — (15 Aug 2024 at 03:00)

Unfortunately, I cannot run your code in MATLAB Online. About two minutes after starting it, it produces:

Error using cat

Requested 4x7x23738800 (5.0GB) array exceeds maximum array size preference (5.0GB). This might cause MATLAB to become

unresponsive.

Error in

ode45 (line 425)

f3d = cat(3,f3d,zeros(neq, 7, chunk, 'like', prototype));

Error in

Alfred_2024_08_15>toODE (line 48)

sol = ode45(@(t,y)batchferment(t,y,b),tspan,y0);

Error in optim.problemdef.fcn2optimexpr

Error in optim.problemdef.fcn2optimexpr

Error in fcn2optimexpr

Error in

Alfred_2024_08_15 (line 56)

myfcn = fcn2optimexpr(@toODE,b,timeex,y0);

Since you will likely be able to run your code (I am not certain I would be able to, for the same reason MATLAB Online cannot), please post back any subsequent errors that your code may throw. It may be possible to correct them without actually running your code.

.

[Alfred]

Okay. Thank you @Star Strider. I will try again.

[Star Strider]

My pleasure!

The problem is likely machine limitation on my end. With my changes, you code works, however it is simply too large for MATLAB Online. Let me know if you have problems.

Getting the correct parameter estimates may not be straightforward with large problems using lsqcurvefit. In that instance, the ga (genetic algoriithm) function is an option if you have tthe Global Optimization Toolbox. I can help with that. The best part about ga is that if there is a problem integrating your differential equations, in most instances, ga just starts over rather than crashing and stopping.

[Alfred]

I have the Global Optimization Tool. You can help me try it.

I am a beginner using MATLAB and need to use it for my project, so I have to learn how to use it.

If you can educate me more, that would be helpful.

Thank you @Star Strider

[Star Strider]

As always, my pleasure!

My revision of your code works. The problem is that it produces matrices that are too large for MATLAB Online.

Does it run successfully on your computer?

EDIT — (26 Aug 2024 at 17:37)

This is not a problem-oriented approach (that may not have existed when I wrote this code), however to use the Genetic Algoriithm (ga) for projects such as these, this approach works —

t=[0.1
0.2
0.4
0.6
0.8
1
1.5
2
3
4
5
6];

c=[0.902	0.06997	0.02463	0.00218
0.8072	0.1353	0.0482	0.008192
0.6757	0.2123	0.0864	0.0289
0.5569	0.2789	0.1063	0.06233
0.4297	0.3292	0.1476	0.09756
0.3774	0.3457	0.1485	0.1255
0.2149	0.3486	0.1821	0.2526
0.141	0.3254	0.194	0.3401
0.04921	0.2445	0.1742	0.5277
0.0178	0.1728	0.1732	0.6323
0.006431	0.1091	0.1137	0.7702
0.002595	0.08301	0.08224	0.835]; 
% theta0=[1;1;1;1;1;1];

% % [theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
ftns = @(theta) norm(c-kinetics(theta,t));
PopSz = 150;
Parms = 10;

optsAns = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10);                                                                  % Options Structure For 'Answers' Problems
% opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);

t0 = clock;
fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)

Start Time: 2024-08-26 17:30:18.0197

[theta,fval,exitflag,output,population,scores] = ga(ftns, Parms, [],[],[],[],zeros(Parms,1),Inf(Parms,1),[],[],optsAns);

ga stopped because the average change in the fitness value is less than options.FunctionTolerance.

t1 = clock;
fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)

Stop Time: 2024-08-26 17:31:25.5311

GA_Time = etime(t1,t0)

GA_Time = 67.5114

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)

Elapsed Time:   6.751140800000000E+01		00:01:07.511

fprintf('Fitness value at convergence = %.4f\nGenerations \t\t\t\t = %d\n\n',fval,output.generations)

Fitness value at convergence = 0.2328
Generations 				 = 1616

fprintf(1,'\tRate Constants:\n')

	Rate Constants:

for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%2d) = %8.5f\n', k1, theta(k1))
end

		Theta( 1) =  0.03251
		Theta( 2) =  1.10289
		Theta( 3) =  4.60035
		Theta( 4) = 10.08034
		Theta( 5) =  1.05261
		Theta( 6) =  0.30051
		Theta( 7) =  1.02803
		Theta( 8) =  0.05707
		Theta( 9) =  0.02482
		Theta(10) =  0.00001

tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);

figure
hd = plot(t, c, 'p');
for k1 = 1:size(c,2)
    CV(k1,:) = hd(k1).Color;
    hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Cfit);
for k1 = 1:size(c,2)
    hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, compose('C_%d',1:size(c,2)), 'Location','N')

function C=kinetics(theta,t)
c0 = theta(7:10);
% c0=[1;0;0;0];
[T,Cv]=ode45(@DifEq,t,c0);
%
    function dC=DifEq(t,c)
    dcdt=zeros(4,1);
    dcdt(1)=-theta(1).*c(1)-theta(2).*c(1);
    dcdt(2)= theta(1).*c(1)+theta(4).*c(3)-theta(3).*c(2)-theta(5).*c(2);
    dcdt(3)= theta(2).*c(1)+theta(3).*c(2)-theta(4).*c(3)+theta(6).*c(4);
    dcdt(4)= theta(5).*c(2)-theta(6).*c(4);
    dC=dcdt;
    end
C=Cv;
end

Use the parameter estimates (‘theta’ values) as initial parameter estimates to lsqcurvefit if they give a good fit to your data with ga. (I have to use ‘optsAns’ to run it here, however using the ‘opts’ option structure will give you more information when you run it on your own computer.)

.

[Alfred]

I have refined some parts, but still errors.

Find the codes that I had run from the same Excel data.

@Star Strider

@Torsten

Kindly help me diagnose the problem. I don't understand the errors.

Thank you.

%Fermentation data
Xdata = readmatrix('Batch1.xlsx','Sheet','sheet1','Range','e2:e13');
Gdata = readmatrix('Batch1.xlsx','Sheet','sheet1','Range','b2:b13');
Bdata = readmatrix('Batch1.xlsx','Sheet','sheet1','Range','c2:c13');
Edata = readmatrix('Batch1.xlsx','Sheet','sheet1','Range','d2:d13');

% Group all data into yariable y
yinv =[Xdata'; Gdata'; Bdata'; Edata'];

%Data for time
timeex = readmatrix('Batch1.xlsx','Sheet','sheet1','Range','a2:a13');

%Set ODE solver
% Initial cell biomass concentration (Initial condition)
y0=[1.04 0.78 0.00 0.00];

%Fermentation time
tspan = (0:2:72);


%Set optimization problem
%Let b = matrix of paramters to be determined
% b= [um ks k1 K1G K1B k2 K2G YXG m k3]
b = optimvar('b',9,"LowerBound",0,"UpperBound",20.75);


%Set functions for ODE solver for solving ODE
function solferment = toODE(b,tspan,y0)
    sol = ode45(@(t,y)batchferment(t,y,b),tspan,y0);
    solferment = deval(sol,tspan);
end

%Convert function for ODE solving to optimization expression 
%To use RtoODE in an objective function, convert the function to an
%optimization expression by using fcn2optimexpr. 
myfcn = fcn2optimexpr(@toODE,b,timeex,y0);

%Express the objective function as the sum of squared differences between
%the ODE solution and the solution with true parameters.
SSE = sum(sum((myfcn-yinv).^2));

%Create an optimization problem with the objective function SSE.
prob = optimproblem ("Description","Fit ODE parameters",'ObjectiveSense','min');

%Objective function (to be minimized)
prob.Objective = SSE;

%Show structure of problem
show(prob)

%Solve Problem 
%To find the best-fitting parameters x, give an initial guess
%x0 for the solver and call solve. 
% Set initial guesses of parameters
%initialGuess.b = [um ks k1 K1G K1B k2 K2G m k3]
initialGuess.b = [0.18 1.0 0.61 0.18 5.85 3.20 16.25 3.02 3.02];

%Solve optimization problem
[sol,optval] = solve(prob,initialGuess);

%Extract the results
%Fitted Parameters
bfinal =sol.b;

%Sum of square Error
SSEfinal = optval;

%Plot the simulated data and experimental data
%Call ODE to solve an equation using Final Fitted Parameters (bfinal)
solysim = ode45(@(t,y)batchferment(t,y,bfinal),tspan,y0);

%Evaluate the simulated results at specified tspan
ysim = deval(solysim,tspan);

%Plot graphs
figure;
plot(timeex,Xdata,'*r',tspan,ysim(1,:),'-b')
legend('exp','sim')
xlabel('Time (h)')
ylabel('X (g/L)')
figure;
plot(timeex,Gdata','*r',tspan,ysim(2,:),'-b')
legend('exp','sim')
xlabel('Time (h)')
ylabel('G (g/L)')
figure;
plot(timeex,Bdata','*r',tspan,ysim(3,:),'-b')
legend('exp','sim')
xlabel('Time (h)')
ylabel('B (g/L)')
figure;
plot(timeex,Edata','*r',tspan,ysim(4,:),'-b')
legend('exp','sim')
xlabel('Time (h)')
ylabel('E (g/L)')

%Equations for Batch 
%y(1) = X = Biomass Concentration (g/l)
%y(2) = G = Glucose Concentration (g/l)
%y(3) = B = Cellobiose Concentration (g/L)
%y(4) = E = Ethanol Concentration (g/l)
function dydt = batchferment(t,y,b)
   % y(1) = X; 
   % y(2) = G;
   % y(3) = B;   
   % y(4) = E;
   X = y(1);
   G = y(2);
   B = y(3);
   E = y(4);

    %%%Growth equations
    %dx/dt = (um*X*G)/(KG+G)
    u = b(1)*y(1)/(b(2)+y(1));
    
    %parameters
  % b(1) = um (1/h);
  % b(2) = ks (g/L);
  um = b(1);
  ks = b(2);

    %%%Cellobiose equations
    C = 88.95-(0.9*y(2))-(0.947*y(3))-(0.9*y(4)/0.511)-(1.137*(y(1)-1.04));
     r1 = (b(3)*C)/(1+(y(2)/b(4))+(y(3)/b(5)));
     r2 = (b(6)*y(3))/(1+(y(2)/b(7)));
  % b(3) = k1;
  % b(4) = K1G;
  % b(5) = K1B;
  % b(6) = k2;
  % b(7) = K2G;
  k1 = b(3);
  K1G = b(4);
  K1B = b(5);
  k2 = b(6);
  k2G = b(7);
  
  %C = Cellulose concentration (g/l)
  %C = C0-(0.9*G)-(0.947*B)-(0.9*E/0.511)-(1.137*(X-X0))
  %r1 = (k1*C)/(1+(G/K1G)+(B/K1B))
  %r2 = (k2*B)/(1+(G/K2G)
 %C = C0-(0.9*G)-(0.947*B)-(0.9*E/0.511)-(1.137*(X-X0))
  %C = 83.65-(0.9*y(2))-(0.947*y(3))-(0.9*y(4)/0.511)-(1.137*(y(1)-1.04))
  
    
    %%% Glucose equations
   rG = ((1/0.11)*(u*y(1)))-(b(8)*y(1));
 
   % b(8) = m;
  m = b(8);

  %rG = ((1/YXG)*(dX/dt))-(m*X)
    
    %%% Ethanol Concentration
    
    % b(9) = k3;
    k3 = b(9); 
    Et = b(9)*(u*y(1))*(1/0.11);  

    %Material balance equations
    dXdt = u*y(1); 
    dGdt = (r2/0.95)-rG;
    dBdt = (r1/0.947)-r2;
    dEdt = Et;
    dydt = [dXdt;dGdt;dBdt;dEdt];
end

The errors are below:

Warning: Failure at t=7.721088e-01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.776357e-15) at time

t.

> In ode45 (line 350)

In Trial>toODE (line 29)

In optim.problemdef.fcn2optimexpr

In optim.problemdef.fcn2optimexpr

In fcn2optimexpr

In Trial (line 36)

Error using deval (line 139)

Attempting to evaluate the solution outside the interval [0.000000e+00, 7.721088e-01] where it is defined.

Error in Trial>toODE (line 30)

solferment = deval(sol,tspan);

Error in optim.problemdef.fcn2optimexpr

Error in optim.problemdef.fcn2optimexpr

Error in fcn2optimexpr

Error in Trial (line 36)

myfcn = fcn2optimexpr(@toODE,b,timeex,y0);

Caused by:

Function evaluation failed while attempting to determine output size. The function might contain an error, or might not be well-defined at the

automatically-chosen point. To specify output size without function evaluation, use 'OutputSize'.

>>

[Alfred]

Hi @Star Strider

https://www.mathworks.com/matlabcentral/answers/2145384-how-to-fit-kinetic-model-to-estimate-kinetic-parameter-from-experimental-data#comment_3246519

Thank you. The code worked on my computer. However, as I am new to MATLAB, I don't fully understand the code you’ve written. Could you kindly explain it to me?

Thank you, and I look forward to your reply.

[Alfred]

Additionally, I noticed that the axis lengths in the output graph are different from what I expected. For example, the maximum time should be 72 hours, and the concentration should be 18 g/L, but the output graph shows different values. Could you explain to me, especially the steps taking throughout?

Thank youv@Star Strider

[Star Strider]

As always, my pleasure!

The essence of that code is that it creates an objective function that ga uses in its fitness function, ‘ftns’. The objective function gets the parameter vector ‘theta’, and independent variable vector (in this instance, time ‘t’) as its only arguments, and uses those to integrate the differential equation system ‘DifEq’ that it then returns to the optimisation function (ga here, it could be others) as the output ‘C’. (In other situations, the objective function might simply be an equation of some sort, however here it is the integrated system of differential equations, since in most instances the system is nonlinear and cannot be solved analytically.) Also, ‘C’ does not need to return all the columns of ‘Cv’ although it does here. If you only had one concentration to fit, you could return only that column of ‘Cv’. The code would still integrate all the equations and use all the parameters.

The optimisation function then compares the results of the integrated differential equations with the data, and adjusts the parameters until it gets a good fit to the data.

The rest of the code simply presents the fitted parameters and plots the results.

The time vector should have specific values, and every element of it should correspond to a specific row of the data matrix you are intending to fit to your system of differential equations. The last time I ran your code I could not get the system to converge. If you are using different code now, please post it and I will see if I can run it. I will then be able to answer your questions about it.

[Alfred]

Thank you @Star Strider

However, I have no knowledge on how to use the Genetic Algoriithm (ga) for projects like this, though this approach works.

I have run a similar code, but I have some errors. Please look at it.

[Alfred]

I have attached the data and the matlab code. When I run the codes, but it seems something is wrong. The graphs look so different that what I expected. Kindly look at it for me.

Thank you @Star Strider

run('Newcode1.mlx')

prob = 
  OptimizationProblem with properties:

       Description: 'Fit ODE parameters'
    ObjectiveSense: 'minimize'
         Variables: [0x0 struct] containing 0 OptimizationVariables
         Objective: [0x0 OptimizationExpression]
       Constraints: [0x0 struct] containing 0 OptimizationConstraints

  No problem defined.

prob = 
  OptimizationProblem with properties:

       Description: 'Fit ODE parameters'
    ObjectiveSense: 'minimize'
         Variables: [1x1 struct] containing 1 OptimizationVariable
         Objective: [1x1 OptimizationExpression]
       Constraints: [0x0 struct] containing 0 OptimizationConstraints

  See problem formulation with show.

  OptimizationProblem : Fit ODE parameters

	Solve for:
       b

	minimize :
       sum(sum((arg1 - extraParams{3}).^2))

       where:

           arg1 = toODE(b, extraParams{1}, extraParams{2});

         extraParams{1}:

          0
          2
          4
          6
          8
         12
         18
         24
         36
         48
         60
         72
         84
         96


       extraParams{2}:

         0.7800         0         0    1.0400


       extraParams{3}:

       Columns 1 through 18

         0.7758         0         0    1.0424    7.4448    1.4314    3.1896    1.0508   11.7505    2.8813    3.4608    1.3110   15.9890    3.0877    3.5559    1.3165   20.0263
         1.3165

       Columns 19 through 20

         4.6822    1.2770

         :
         :




	variable bounds:
       0 <= b(1)  <= 25
       0 <= b(2)  <= 25
       0 <= b(3)  <= 25
       0 <= b(4)  <= 25
       0 <= b(5)  <= 25
       0 <= b(6)  <= 25
       0 <= b(7)  <= 25
       0 <= b(8)  <= 25
       0 <= b(9)  <= 25
       0 <= b(10) <= 25


Solving problem using lsqnonlin.

Solver stopped prematurely.

lsqnonlin stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 1.000000e+03.

sol = struct with fields:
    b: [10x1 double]

optval = 1.1132e+60

bfinal = 10×1

0.2061 1.0014 0.6077 0.1752 5.7602 3.1736 16.2326 0.1350 3.3820 3.0119

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

SSEfinal = 1.1132e+60

solYsim = struct with fields:
     solver: 'ode45'
    extdata: [1x1 struct]
          x: [0 2.3401e-05 1.4041e-04 7.2543e-04 0.0037 0.0183 0.0914 0.3033 0.5914 0.9210 0.9724 1.0066 1.0176 1.0218 1.0260 1.0282 1.0285 1.0286 1.0286 1.0287 ... ] (1x102 double)
          y: [4x102 double]
      stats: [1x1 struct]
      idata: [1x1 struct]

Ysim = 4×100

1.0e+29 * 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0001 -0.0004 -0.0015 -0.0046 -0.0126 -0.0302 -0.0643 -0.1226 0 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0004 0.0014 0.0044 0.0120 0.0287 0.0610 0.1165 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Alfred]

Any feedback? @Star Strider @Torsten

Thank you.

[Torsten]

It would help if you included the .m-file instead of the .mlx-file.

[Alfred]

ok. Here is the .m file. @Torsten

[Star Strider]

I ran this in MATLAB Online, and it essentially works, in the sense that it runs without any code errors. (I changed the code to correspond to my usual format., since I understand it and I know that it works.)

The problem is that it throws this process error:

Warning: Failure at t=1.113343e+01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.842171e-14) at time t.

The integration does not complete at that point, throwing the second error:

Arrays have incompatible sizes for this operation.

That refers to the comparison of the arrays in the ‘ftns’ function.

I am not certain what the problem is, and since I do not know what equations you are integrating, I cannot check them for codiing errors. (That is always the first thing I do when my own code fails.) At some point (this varies so it is not always at the same time), the integrated values result in a 0/0 or Inf/Inf calculation, and that produces a NaN value in the integrated output. The ga iteration then stops, since the data matrix and the matrix returned by the ‘toODE’ function then have incompatible sizes.

Please check to be certain that the differential equations are coded correctly. My code should then produce optimised parameters.

%Experimantal Biomass data
Sexp = readmatrix('Data.xlsx','Sheet','2','Range','d5:d18');
Bexp = readmatrix('Data.xlsx','Sheet','2','Range','e5:e18');
Pexp = readmatrix('Data.xlsx','Sheet','2','Range','f5:f18');
Xexp = readmatrix('Data.xlsx','Sheet','2','Range','g5:g18');

% Group all data into yariable Y
Yexp = [Sexp Bexp Pexp Xexp];


%Exprimantal time data 
timeexp = readmatrix('Data.xlsx','Sheet','2','Range','c5:c18');

ftns = @(b) norm(Yexp - toODE(b,timeexp));
PopSz = 150;
Parms = 14;

% optsAns = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10);                                                                  % Options Structure For 'Answers' Problems
opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-5, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);

t0 = clock;
fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)
[theta,fval,exitflag,output,population,scores] = 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('Fitness value at convergence = %.4f\nGenerations \t\t\t\t = %d\n\n',fval,output.generations)

fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%2d) = %8.5f\n', k1, theta(k1))
end

tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);

figure
hd = plot(t, c, 'p');
for k1 = 1:size(c,2)
    CV(k1,:) = hd(k1).Color;
    hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Cfit);
for k1 = 1:size(c,2)
    hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, compose('C_%d',1:size(c,2)), 'Location','N')

function solferment = toODE(b,tspan)
% fprintf(['\nb1 = ',repmat('%f ', 1, 14) '\n'],b)
Y0 = b(end-3:end);
% function solferment = toODE(b,tspan,Y0)
%call ODE45
[t,Cv] = ode15s(@(t,Y)fermenterODEs(t,Y,b),tspan,Y0);
%Evaluate solution using tspan
% % solferment = deval(sol,tspan);
end
%b = [um ks Pmax nx qpm ksp Ppmax np Yxs Yps] 
%b(1) = um 
%b(2) = ks 
%b(3) = k1 
%b(4) = K1G 
%b(5) = K1B 
%b(6) = k2 
%b(7) = K2G 
%b(8) = YXG 
%b(9) = m 
%b(10) = k3     

function dYdt = fermenterODEs(t,Y,b)
% function dYdt = fermenterODEs(t,Y,b)
% fprintf(['\nb2 = ',repmat('%f ', 1, 10) '\n'],b)
    S = Y(1); % Sustrate, S
    B = Y(2); % Cellobiose, B
    P = Y(3); % Product, S
    X = Y(4); % Biomass, x


    %Growth Kinetic equations
    u =(b(1)*Y(1)/(b(2)+Y(1)));

    %dX/dt = u*Y(1);


%%%% Glucose concentration

%dG/dt = (r2/0.95)-rG;

    rG = ((u*Y(4))/b(8))+(b(9)*Y(4));


    %Cellobiose equations
     %C0 = 73.658 g/L cellulose concentration
     %X0 = 0.18 g/L biomass concentration
  
     %%%% Cellobiose Concentration

  C0 = 73.66;
  X0 = 1.04;
  
  C = C0-(0.9*Y(1))-(0.947*Y(2))-(0.9*Y(3)/0.511)-(1.137*(Y(4)-X0));
  %C = 73.658-(0.9*y(3))-(0.947*y(4))-(0.9*y(2)/0.511)-(1.137*(y(1)-0.18))
  %C = b(20)-(0.9*y(3))-(0.947*y(4))-(0.9*y(2)/0.511)-(1.137*(y(1)-0.18))
      r1 = (b(3)*C)/(1+(Y(2)/b(5))+(Y(1)/b(4)));
     r2 = (b(6)*Y(2))/(1+(Y(1)/b(7)));

   %dB/dt = (r1/0.947)-r2; 

%%%% Ethanol Concentration
  %dE/dt = (b(10)*b(1)*y(4)*y(1))/(b(2)+y(1))*b(8);

  

    %Material balance equations
    dSdt = u*Y(1);
    dBdt = (r2/0.95)-rG;
    dPdt = (r1/0.947)-r2; 
    dXdt = (b(10)*b(1)*Y(4)*Y(1))/(b(2)+Y(1))*b(8);
    
    dYdt = [dSdt;dBdt;dPdt;dXdt];
end

II will keep my MATLAB Online code up so I can run it again with the correct differential equations. I plotted the data, and while this may be a challenge to optimise (it may be neceessary to run my code several times before it finds an optimal set of parameters), that should be possible.

NOTE — As is my usual practice, I included the initial condition vector ‘Y0’ separately as the last four (in this instance) elements in the parameter vector, so the parameter vector has a total of 14 elements. I always estimate the initial conditions as well as the desired parameters because that produces a more reliable result. The initial values in the data are themselves measured quantities, and so subject to error.

.

[Torsten]

Before starting an optimization, you should at least make sure that you get reasonable curves for your initial parameter vector b. Now look below for the curves produced for your initial b-vector. An optimizer has no chance to succeed with these initial simulated curves for S,B,P and X.

%% Fitting ODE Parameters (Batch Fermentation)
% Load experimental data

%Experimantal Biomass data
Sexp = readmatrix('Data.xlsx','Sheet','2','Range','d5:d18');
Bexp = readmatrix('Data.xlsx','Sheet','2','Range','e5:e18');
Pexp = readmatrix('Data.xlsx','Sheet','2','Range','f5:f18');
Xexp = readmatrix('Data.xlsx','Sheet','2','Range','g5:g18');

% Group all data into yariable Y
Yexp =[Sexp'; Bexp'; Pexp'; Xexp'];

%Exprimantal time data 
timeexp = readmatrix('Data.xlsx','Sheet','2','Range','c5:c18'); 
% Set  ODE solver
% Define initial condition for ODE solver

%initial biomass concentration
Y0 = [0.78 0.00 00.00 1.04];
%% 
% Define time range for ODE solver
b = [0.18 1.00 0.61 0.18 5.85 3.20 16.25 0.11 3.40 3.02];
%
tspan = timeexp(1):timeexp(end);
%Call ODE to solve an equation using Final Fitted Parameters (bfinal)
[Tsim,Ysim] = ode15s(@(t,Y)fermenterODEs(t,Y,b),tspan,Y0);
Ysim = Ysim.';
% Plot graphs
plot(timeexp,Xexp,'*r',tspan,Ysim(1,:),'-b')
xlabel('Time (h)')
ylabel('S (g/L)')

plot(timeexp,Pexp','*r',tspan,Ysim(2,:),'-b')
xlabel('Time (h)')
ylabel('C (g/L)')

plot(timeexp,Sexp','*r',tspan,Ysim(3,:),'-b')
xlabel('Time (h)')
ylabel('P (g/L)')

plot(timeexp,Sexp','*r',tspan,Ysim(3,:),'-b')
xlabel('Time (h)')
ylabel('X (g/L)')

% Set functions for ODE solver
% This function is for solving ODE


% This function is for writing  equation

%b = [um ks Pmax nx qpm ksp Ppmax np Yxs Yps] 
%b(1) = um 
%b(2) = ks 
%b(3) = k1 
%b(4) = K1G 
%b(5) = K1B 
%b(6) = k2 
%b(7) = K2G 
%b(8) = YXG 
%b(9) = m 
%b(10) = k3     

function dYdt = fermenterODEs(t,Y,b)

    S = Y(1); % Sustrate, S
    B = Y(2); % Cellobiose, B
    P = Y(3); % Product, S
    X = Y(4); % Biomass, x


    %Growth Kinetic equations
    u =(b(1)*Y(1)/(b(2)+Y(1)));

    %dX/dt = u*Y(1);


%%%% Glucose concentration

%dG/dt = (r2/0.95)-rG;

    rG = ((u*Y(4))/b(8))+(b(9)*Y(4));


    %Cellobiose equations
     %C0 = 73.658 g/L cellulose concentration
     %X0 = 0.18 g/L biomass concentration
  
     %%%% Cellobiose Concentration

  C0 = 73.66;
  X0 = 1.04;
  
  C = C0-(0.9*Y(1))-(0.947*Y(2))-(0.9*Y(3)/0.511)-(1.137*(Y(4)-X0));
  %C = 73.658-(0.9*y(3))-(0.947*y(4))-(0.9*y(2)/0.511)-(1.137*(y(1)-0.18))
  %C = b(20)-(0.9*y(3))-(0.947*y(4))-(0.9*y(2)/0.511)-(1.137*(y(1)-0.18))
      r1 = (b(3)*C)/(1+(Y(2)/b(5))+(Y(1)/b(4)));
     r2 = (b(6)*Y(2))/(1+(Y(1)/b(7)));

   %dB/dt = (r1/0.947)-r2; 

%%%% Ethanol Concentration
  %dE/dt = (b(10)*b(1)*y(4)*y(1))/(b(2)+y(1))*b(8);

  

    %Material balance equations
    dSdt = u*Y(1);
    dBdt = (r2/0.95)-rG;
    dPdt = (r1/0.947)-r2; 
    dXdt = (b(10)*b(1)*Y(4)*Y(1))/(b(2)+Y(1))*b(8);
    
    dYdt = [dSdt;dBdt;dPdt;dXdt];
end

[Alfred]

Thank you @Star Strider

Based on the comment you made: https://www.mathworks.com/matlabcentral/answers/2145384-how-to-fit-kinetic-model-to-estimate-kinetic-parameter-from-experimental-data#comment_3295390

I think you should see the equations. I am certain that the differential equations are coded correctly.

You can have a look at the equation in the attached file and the paper from with I got the equations.

What do you think, to try to fit the data first and the optimize the parameters or do them together? Can help me fit the data and see?

Thank you very much.

[Alfred]

Thank you @Torsten

Based on your comment: https://www.mathworks.com/matlabcentral/answers/2145384-how-to-fit-kinetic-model-to-estimate-kinetic-parameter-from-experimental-data#comment_3295410

What can I do? How can I get reasonable curves for your initial parameter vector b?

I am a novice novice in Matlab, but I need to use it in my work as soon as possible, that's why I am striving.

Thank you.

You can find the equations in the files below with a paper from with I got the equations.

[Star Strider]

My pleasure!

In the provided document, there are five diferential equations. In your code, what variables correspond to ‘C’, ‘B’, ‘G’, ‘X’, and ‘E’?

That information will make it easier for me to follow this.

[Alfred]

Thank you for your rapid reply @Star Strider.

In the file named "Kinetic modelling_Alfred,"

‘C’ = Cellulose

‘B’ = Cellobiose

‘G’ = Glucose

‘X’ = Cell biomass

‘E’ = Ethanol

[Star Strider]

Thank you.

However in your latest code, the variiables are:

S = Y(1); % Sustrate, S

B = Y(2); % Cellobiose, B

P = Y(3); % Product, S

X = Y(4); % Biomass, x

and there are four differential equations:

%Material balance equations

dSdt = u*Y(1);

dBdt = (r2/0.95)-rG;

dPdt = (r1/0.947)-r2;

dXdt = (b(10)*b(1)*Y(4)*Y(1))/(b(2)+Y(1))*b(8);

instead of the five presented in the PDF document. I need to understand this.

(My apologies for the delay. Out to run errands today.)

.

[Torsten]

What can I do? How can I get reasonable curves for your initial parameter vector b?

The question is: How can I get an initial parameter vector b to get reasonable curves at the start ?

If your model is published in a paper, you should start with the parameters used therein as initial guesses. Although the produced curves might be far off, they will at least be reasonable without results in the order of 1e33.

[Alfred]

@Torsten

Yes, that is what I used. I used the parameters published as initial guesses. Yet, getting such curves.

[Star Strider]

This produces output that matches the input matrix, and produces reasonable reaults, however I am having problems getting it to converge.

I re-derived the differenttial equations using your p[rovided PDF and the Symbolic Math Toolbox witth this code:

clear all

close all

syms C(t) B(t) G(t) X(t) E(t) YXG k1prime k1 k2 k3 KG K1B K1G K2G m um T Y

r1 = k1prime * C / (1 + B/K1B + G/K1G)

r2 = k2 * B / (1 + G/K2G)

rx = um * X * G / (K1G + G)

rG = rx/YXG + m * X

Eq1 = diff(C) == -r1

Eq2 = diff(B) == r1/0.947 - r2

Eq3 = diff(G) == r2/0.95 - rG

Eq4 = diff(X) == um * X * G / (KG + G)

Eq5 = diff(E) == k3 * um * X * G / ((KG + G) * YXG)

[VF,Subs] = odeToVectorField(Eq1, Eq2, Eq3, Eq4, Eq5)

DEfcn = matlabFunction(VF, 'Vars',{t,Y,K1B,K1G,K2G,KG,YXG,k2,k3,k1prime,m,um})

t = 1:10;

[t,Y] = ode45(@(t,Y)DEfcn(t,Y,rand,rand,rand,rand,rand,rand,rand,rand,rand,rand),t,rand(1,5))

figure

plot(t,Y)

grid

I invite you to check that to see if there are any errors, correct them as necessary, then run it. (I did not see any errors, however I was working from your PDF, so the code is as good as those equations.) I then used ‘DEfcn’ in the code to fit the data. (II adapted my original ‘kinetics’ code for this problem, since I understand it.) The ‘Subs’ comments are returned by the odeToVectorField function, and correspond to the order of the differential equations in ‘DEfcn’.

I integrated the code using random parameters to be certain that it produces results resembling the original data, and it seems to provide a reasonable result. I just have problems fitting it.

So I invite you to check it, make any necessary corrections, and post back with results.

The ‘Arrays have incompatible sizes for this operation’ error means that the integratiion encountered a singluarity and did not complete, so the comparison inn the ‘ftns’ function fails. If it does not encounter the singularity, the matrices will be the same size, and the code will continue.

%Fermentation data

Xdata = readmatrix('Batch1.xlsx','Sheet','sheet1','Range','e2:e13');
Gdata = readmatrix('Batch1.xlsx','Sheet','sheet1','Range','b2:b13');
Bdata = readmatrix('Batch1.xlsx','Sheet','sheet1','Range','c2:c13');
Edata = readmatrix('Batch1.xlsx','Sheet','sheet1','Range','d2:d13');


% Group all data into yariable y
yinv =[Xdata Gdata Bdata Edata];
c = yinv;

%Data for time
timeex = readmatrix('Batch1.xlsx','Sheet','sheet1','Range','a2:a13');
t = timeex;



% % [theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
ftns = @(theta) norm(c-kinetics(theta,t));
PopSz = 50;
Parms = 15;

optsAns = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10);                                                                  % Options Structure For 'Answers' Problems
% opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);

t0 = clock;
fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)

Start Time: 2024-11-09 21:27:54.7683

[theta,fval,exitflag,output,population,scores] = ga(ftns, Parms, [],[],[],[],zeros(Parms,1),Inf(Parms,1),[],[],optsAns);

Warning: Failure at t=3.009454e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (7.105427e-15) at time t.

Arrays have incompatible sizes for this operation.

Error in solution>@(theta)norm(c-kinetics(theta,t)) (line 20)
ftns = @(theta) norm(c-kinetics(theta,t));

Error in createAnonymousFcn>@(x)fcn(x,FcnArgs{:}) (line 11)
fcn_handle = @(x) fcn(x,FcnArgs{:});

Error in makeState (line 58)
firstMemberScore = FitnessFcn(state.Population(initScoreProvided+1,:));

Error in galincon (line 24)
state = makeState(GenomeLength,FitnessFcn,Iterate,output.problemtype,options);

Error in ga (line 420)
[x,fval,exitFlag,output,population,scores] = galincon(FitnessFcn,nvars, ...

Caused by:
Failure in initial user-supplied fitness function evaluation. GA cannot continue.

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('Fitness value at convergence = %.4f\nGenerations \t\t\t\t = %d\n\n',fval,output.generations)

fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%2d) = %8.5f\n', k1, theta(k1))
end

tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);

figure
hd = plot(t, c, 'p');
for k1 = 1:size(c,2)
    CV(k1,:) = hd(k1).Color;
    hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Cfit);
for k1 = 1:size(c,2)
    hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, compose('C_%d',1:size(c,2)), 'Location','N')

t = timeex;

function C=kinetics(theta,t)
c0 = theta(11:15);
% [um,KG,k1prime,K1G,K1B,k2,K2G,YXG,m,k3] = num2cell(theta(1:10))
um = theta(1);
KG = theta(2);
k1prime = theta(3);
K1G = theta(4);
K1B = theta(5);
k2 = theta(6);
K2G = theta(7);
YXG = theta(8);
m = theta(9);
k3 = theta(10);

DEfcn = @(t,Y,K1B,K1G,K2G,KG,YXG,k2,k3,k1prime,m,um)[-(k2.*Y(1))./(Y(3)./K2G+1.0)+(k1prime.*Y(2).*(1.0e+3./9.47e+2))./(Y(1)./K1B+Y(3)./K1G+1.0);-(k1prime.*Y(2))./(Y(1)./K1B+Y(3)./K1G+1.0);-m.*Y(4)+(k2.*Y(1).*(2.0e+1./1.9e+1))./(Y(3)./K2G+1.0)-(um.*Y(3).*Y(4))./(YXG.*(K1G+Y(3)));(um.*Y(3).*Y(4))./(KG+Y(3));(k3.*um.*Y(3).*Y(4))./(YXG.*(KG+Y(3)))];

% Y = rand(5,1);
% 
% QQ = DEfcn(t,Y,K1B,K1G,K2G,KG,YXG,k2,k3,k1prime,m,um)

[t,Cv] = ode15s(@(t,Y)DEfcn(t,Y,K1B,K1G,K2G,KG,YXG,k2,k3,k1prime,m,um),t,c0);

C = Cv(:,[4,3,1,5]);

% Subs =
% 
% B
% C
% G
% X
% E

end

% return
% 
% % Subs =
% % 
% % B
% % C
% % G
% % X
% % E
% % c0 = theta(7:10);
% % % c0=[1;0;0;0];
% % [T,Cv]=ode45(@DifEq,t,c0);
% % %
% %     function dC=DifEq(t,c)
% %     dcdt=zeros(4,1);
% %     dcdt(1)=-theta(1).*c(1)-theta(2).*c(1);
% %     dcdt(2)= theta(1).*c(1)+theta(4).*c(3)-theta(3).*c(2)-theta(5).*c(2);
% %     dcdt(3)= theta(2).*c(1)+theta(3).*c(2)-theta(4).*c(3)+theta(6).*c(4);
% %     dcdt(4)= theta(5).*c(2)-theta(6).*c(4);
% %     dC=dcdt;
% %     end
% % C=Cv;
% 
% 
% return
% 
% %Set ODE solver
% % Initial cell biomass concentration (Initial condition)
% y0=[1.04 0.78 0.00 0.00];
% 
% %Fermentation time
% tspan = linspace(0,72);
% 
% %Set optimization problem
% %Let b = matrix of paramters to be determined
% % b= [um ks k1 K1G K1B k2 K2G YXG m k3]
% b = optimvar('b',10,"LowerBound",0,"UpperBound",72);
% 
% %Set functions for ODE solver for solving ODE
% function solferment = toODE(b,tspan,y0)
%     sol = ode45(@(t,y)batchferment(t,y,b),tspan,y0);
%     solferment = deval(sol,tspan);
% end
% 
% %Convert function for ODE solving to optimization expression 
% 
% %To use RtoODE in an objective function, convert the function to an
% %optimization expression by using fcn2optimexpr. 
% myfcn = fcn2optimexpr(@toODE,b,timeex,y0);
% 
% %Express the objective function as the sum of squared differences between
% %the ODE solution and the solution with true parameters.
% SSE = sum(sum((myfcn-yinv).^2));
% 
% %Create an optimization problem with the objective function SSE.
% prob = optimproblem ("Description","Fit ODE parameters",'ObjectiveSense','min');
% 
% %Objective function (to be minimized)
% prob.Objective = SSE;
% 
% %Show structure of problem
% show(prob)
% 
% %Solve Problem 
% %To find the best-fitting parameters x, give an initial guess
% %x0 for the solver and call solve. 
% % Set initial guesses of parameters
% initialGuess.b = [0.18 1.0 0.61 0.18 5.85 3.20 16.25 0.11 3.40 3.02];
% 
% %Solve optimization problem
% [sol,optval] = solve(prob,initialGuess);
% 
% %Extract the results
% %Fitted Parameters
% bfinal =sol.b;
% 
% %Sum of square Error
% SSEfinal = optval;
% 
% 
% %Plot the simulated data and experimental data
% %Call ODE to solve an equation using Final Fitted Parameters (bfinal)
% solysim = ode45(@(t,y)batchferment(t,y,bfinal),tspan,y0);
% %Evaluate the simulated results at specified tspan
% ysim = deval(solysim,tspan);
% 
% %Plot graphs
% figure;
% plot(timeex,Xdata,'*r',tspan,ysim(1,:),'-b')
% legend('exp','sim')
% xlabel('Time (h)')
% ylabel('X (g/L)')
% 
% figure;
% plot(timeex,Gdata','*r',tspan,ysim(2,:),'-b')
% legend('exp','sim')
% xlabel('Time (h)')
% ylabel('G (g/L)')
% 
% figure;
% plot(timeex,Bdata','*r',tspan,ysim(3,:),'-b')
% legend('exp','sim')
% xlabel('Time (h)')
% ylabel('B (g/L)')
% 
% figure;
% plot(timeex,Edata','*r',tspan,ysim(4,:),'-b')
% legend('exp','sim')
% xlabel('Time (h)')
% ylabel('E (g/L)')
% 
% 
% %Equations for Batch 
% %y(1) = X = Biomass Concentration (g/l)
% %y(2) = G = Glucose Concentration (g/l)
% %y(3) = B = Cellobiose Concentration (g/L)
% %y(4) = E = Ethanol Concentration (g/l)
% 
% function dydt = batchferment(t,y,b)
% 
%    % y(1) = X; 
%    % y(2) = G;
%    % y(3) = B;   
%    % y(4) = E;
%    X = y(1);
%    G = y(2);
%    B = y(3);
%    E = y(4);
% 
%     %%%Growth equations
%     %dx/dt = (um*X*G)/(KG+G)
%     u = (b(1)*y(1)*y(2))/(b(2)+y(2));
%     %u = (b(1)*X*G)/(b(2)+G);
%     %parameters
%   % b(1) = um (1/h);
%   % b(2) = ks (g/L);
%   um = b(1);
%   ks = b(2);
% 
% 
%     %%%Cellobiose equations
%     C = b(1)-(0.9*y(2))-(0.947*y(3))-(0.9*y(4)/0.511)-(1.137*(y(1)-1.04));
%      r1 = (b(3)*C)/(1+(y(2)/b(4))+(y(3)/b(5)));
%      r2 = (b(6)*y(3))/(1+(y(2)/b(7)));
%   % b(3) = k1;
%   % b(4) = K1G;
%   % b(5) = K1B;
%   % b(6) = k2;
%   % b(7) = K2G;
%   k1 = b(3);
%   K1G = b(4);
%   K1B = b(5);
%   k2 = b(6);
%   k2G = b(7);
% 
%   %C = Cellulose concentration (g/l)
%   %C = C0-(0.9*G)-(0.947*B)-(0.9*E/0.511)-(1.137*(X-X0))
%   %r1 = (k1*C)/(1+(G/K1G)+(B/K1B))
%   %r2 = (k2*B)/(1+(G/K2G)
% 
%  %C = C0-(0.9*G)-(0.947*B)-(0.9*E/0.511)-(1.137*(X-X0))
%   %C = 83.65-(0.9*y(2))-(0.947*y(3))-(0.9*y(4)/0.511)-(1.137*(y(1)-1.04))
% 
% 
%     %%% Glucose equations
%     % rG = ((1/b(8))*(dx/dt))-(b(9)*y(1));
%   % b(8) = YXG;
%   % b(9) = m;
%   YXG = b(8)
%   m = b(9);
%   %rG = ((1/YXG)*(dX/dt))-(m*X)
% 
% 
%     %%% Ethanol Concentration
% 
%     % b(10) = k3;
%     k3 = b(10); 
%     Et = k3*(u*X)*(1/YXG);  
% 
%     %Material balance equations
%     dXdt = u*X;
%     rG = ((1/b(8))*(dXdt))-(b(9)*y(1));
%     dGdt = (r2/0.95)-rG;
%     dBdt = (r1/0.947)-r2;
%     dEdt = Et;
%     dydt = [dXdt;dGdt;dBdt;dEdt];
% end

.

[Alfred]

@Star Strider

I am lost with the codes.

I think I should explain a bit more about my situation, I think this will help.

So, I carried out a fermentation experiment and the result is presented in the Excel file named "Batch1" as attached. The equation I am to use is already published in the paper attached.

What I want to do is to fit the data and try to optimize the parameters.

I finding it difficult to write the codes.

However, the solution you had is very good:https://www.mathworks.com/matlabcentral/answers/2145384-how-to-fit-kinetic-model-to-estimate-kinetic-parameter-from-experimental-data#comment_3246519

The only challenge I have with it is that the range of the x and y axis are different from my data.

@Star Strider

@Torsten

Any assistance would be greatly appreciated, as always.

Thank you.

[Star Strider]

Having the paper helps.

I wiill look through it (that will take a while) and make appropriate suggestions.

[Alfred]

Okay. Thank you @Star Strider. I will be looking out for it.

[Star Strider]

That has to be one of the most confusing papers I have ever read. Part of that may be that this is not an area of my expertise, however it is difficult for me to follow it.

As best I can determine, the four differential equations are for: , , , and , where B is cellbiose, G is glucose, X is cell concentration, and E is ethanol concentration. These are all described in section 6, given by equations . Your equations for them do not (to me) appear to be the same as those described in the paper.

Since I do not understand your data with respect to those equations, the first suggestion I have is to check your equations against those four equations to be certain that they are the same, and that they are appropriate to your data (with respect to time units, concentration units, and mass units, and if not, convert the units as necessary to be consistent and compatible with those in the paper).

The equations in the paper and the equations in your code do not appear to me to be the same, however I may be missing somethiing. Please check that, and make any necessary corrections. Then, run your code again and see if it converges on an appropriate solution.

.

[Alfred]

Thank you very much for your suggestions. I will check them out and provide feedback. @Star Strider

Wishing you a Happy New Year filled with greater heights in all your endeavors. @Star Strider and @Torsten

[Star Strider]

As always, my (our) pleasure!

[Alfred]

Hello @Star Strider

Please look at the codes for me. I am having some errors. I have made all the correction as you've recommended. The equations are now the same as in the paper. However, I am getting some errors.

In works such as mine, Lsqcurvefit and Lsqnonlin are mostly used. However, I don't know how to use them.

Thank you.

[Torsten]

Could you include the .m file, not the .mlx file ?

[Alfred]

Okay, sure.

Here!

[Torsten]

dBdt = (b(8)*C)/0.947*(1+(S/b(10))+(B/b(11)))*(b(9)/1+b(9)*(12)*t)-(b(1)*B/(1+(S/b(2))));

I'm quite sure that several parentheses have been set incorrectly in the definition of the time derivatives, e.g. the bold expression above.

[Alfred]

@Torsten

It is correct. For example b(9) is representative of a parameter.

Also, the number of bracket in each side is equal.

[Star Strider]

When I ran this just now, it is having array size problems. I will explore running this on my offline computer later.

%% Fitting ODE Parameters (Batch Fermentation)
% Load experimental data

%Experimantal Biomass data
Xexp = readmatrix('exodefit.xlsx','Sheet','2','Range','d5:d18');
Pexp = readmatrix('exodefit.xlsx','Sheet','2','Range','e5:e18');
Sexp = readmatrix('exodefit.xlsx','Sheet','2','Range','f5:f18');
Bexp = readmatrix('exodefit.xlsx','Sheet','2','Range','g5:g18');

% Group all data into yariable Y
Yexp =[Xexp'; Pexp'; Sexp'; Bexp'];

%Exprimantal time data 
timeexp = readmatrix('exodefit.xlsx','Sheet','2','Range','c5:c18');
% Set  ODE solver
% Define initial condition for ODE solver

%initial biomass concentration
Y0 = [1.04 0.00 0.78 0.00];
%% 
% Define time range for ODE solver

%Fermentation time
tspan = linspace(0,96);
% Set optimization problem
% Declare b = matrix of paramters to be determined

% b= [k2 k2G kG um m k3 YXG k1 e k1G k1B k4] 

b = optimvar('b',12,"LowerBound",0,"UpperBound",100); 
%% 
% Convert function for ODE solving to optimization expression 

%myfcn = solution from ODE solving
myfcn = fcn2optimexpr(@toODE,b,timeexp,Y0)

Error using cat
Requested 4x7x23738800 (5.0GB) array exceeds maximum array size preference (5.0GB). This might cause MATLAB to become unresponsive.

Error in ode45 (line 421)
f3d = cat(3,f3d,zeros(neq, 7, chunk, 'like', prototype));

Error in solution>toODE (line 109)
sol = ode45(@(t,Y)fermenterODEs(t,Y,b),tspan,Y0);

Error in optim.problemdef.fcn2optimexpr

Error in optim.problemdef.fcn2optimexpr

Error in fcn2optimexpr

Caused by:
Function evaluation failed while attempting to determine output size. The function might contain an error, or might not be well-defined at the automatically-chosen point. To specify output size without function evaluation, use 'OutputSize'.

%% 
% Set sum of squared error as an objective function to be minimized

SSE = sum(sum((myfcn-Yexp).^2));
%% 
% Set optimization problem

prob = optimproblem ("Description","Fit ODE parameters",'ObjectiveSense','min')
%% 
% Objective function (to be minimized)

prob.Objective = SSE
%% 
% Show structure of problem

show(prob)
%% 
% Set initial guesses of parameters

%b = [k2 k2G kG um m k3 YXG]
initialGuess.b = [3.20 16.25 3.40 0.18 0.64 3.02 0.52 0.61 1.00 0.18 5.85 0.10];
%b(1) = k2
%b(2) = k2G
%b(3) = kG
%b(4) = um
%b(5) = m
%b(6) = k3
%b(7) = YXG
%b(8) = k1
%b(9) = e
%b(10) = k1G
%b(11) = k1B
%b(12) = k4
%% 
% Solve  optimization problem

[sol,optval] = solve(prob,initialGuess)
%% 
% Extract the results

%Fitted Parameters
bfinal =sol.b
%Sum of square Error
SSEfinal = optval
% Plot the simulated data and experimental data

%Call ODE to solve an equation using Final Fitted Parameters (bfinal)
solYsim = ode45(@(t,Y)fermenterODEs(t,Y,bfinal),tspan,Y0)
%Evaluate the simulated results at specified tspan
Ysim = deval(solYsim,tspan)
%% 
% Plot graphs

plot(timeexp,Xexp,'*r',tspan,Ysim(1,:),'-b')
xlabel('Time (h)')
ylabel('X (g/L)')

plot(timeexp,Pexp','*r',tspan,Ysim(2,:),'-b')
xlabel('Time (h)')
ylabel('P (g/L)')

plot(timeexp,Sexp','*r',tspan,Ysim(3,:),'-b')
xlabel('Time (h)')
ylabel('S (g/L)')

plot(timeexp,Bexp','*r',tspan,Ysim(4,:),'-b')
xlabel('Time (h)')
ylabel('B (g/L)')
% Set functions for ODE solver
% This function is for solving ODE

function solferment = toODE(b,tspan,Y0)
%call ODE45
sol = ode45(@(t,Y)fermenterODEs(t,Y,b),tspan,Y0);
%Evaluate solution using tspan
solferment = deval(sol,tspan);
end
%% 
% This function is for writing  equation

%b = [k2 k2G kG um m k3 YXG]
%b(1) = k2
%b(2) = k2G
%b(3) = kG
%b(4) = um
%b(5) = m
%b(6) = k3
%b(7) = YXG
%b(8) = k1
%b(9) = e
%b(10) = k1G
%b(11) = k1B
%b(12) = k4

function dYdt = fermenterODEs(t,Y,b)

    X = Y(1); % Biomass, x
    P = Y(2); % Product, P
    S = Y(3); % Sustrate, S
    B = Y(4); % Cellobiose, B

    %Growth Kinetic equations
    

    %Product formation kineticss
%b(1) = k2
%b(2) = k2G
%b(3) = kG
%b(4) = um
%b(5) = m
%b(6) = k3
%b(7) = YXG
%b(8) = k1
%b(9) = e
%b(10) = k1G
%b(11) = k1B
%b(12) = k4

    %Material balance equations
    C=40-(0.9*S)-(0.947*B)-(0.9*P/0.511)-1.137*(X-1.04);
    dXdt = (b(4)*X*S)/(b(3)+S);
    dPdt = (b(6)*b(4)*X*S)/(b(3)+S)*b(7);
    dSdt = ((b(1)*B)/(0.95*(1+(S/b(2)))))-((b(4)*X*S)/(b(3)+S)*b(7))-(b(5)*X);
    dBdt = (b(8)*C)/(0.947*(1+(S/b(10))+(B/b(11))))*(b(9)/(1+b(9)*(12)*t))-((b(1)*B)/(1+(S/b(2)))); 
    dYdt = [dXdt;dPdt;dSdt;dBdt];
end
%(b(9)/(1+b(9)*(12)*t))

.

EDIT — (14 Jan 2025 at 17:12)

New script, same problem. I’ll follow up later.

.

[Alfred]

Okay. I will be looking forward to it.

Thank you @Star Strider

[Star Strider]

As always, my pleasure!

I have to be away for a while. I will get back to it when II return.

[Torsten]

@Alfred

It is correct. For example b(9) is representative of a parameter.

Also, the number of bracket in each side is equal.

Shouldn't the expression read

(b(9)/(1+b(9)*(12)*t))

instead of

(b(9)/1+b(9)*(12)*t)

?

If not, why do you divide b(9) by 1 ? b(9)/1 = b(9).

[Alfred]

@Torsten

You are right. I just noticed. Thank you.

[Alfred]

@Star Strider @Torsten

Please find the corrected file.

Thank you.

[Torsten]

You have to start with initial values for b that at least keep the concentrations positive. That's not the case for your choice for b.

Did you manage to use the data from the article and reproduce the respective curves ?

%initial biomass concentration
Y0 = [1.04 0.00 0.78 0.00];
%% 
% Define time range for ODE solver
%Fermentation time
b = [3.20 16.25 3.40 0.18 0.64 3.02 0.52 0.61 1.00 0.18 5.85 0.10];
[T,Y] = ode45(@(t,Y)fermenterODEs(t,Y,b),[0 96],Y0);
plot(T,Y)

function dYdt = fermenterODEs(t,Y,b)
    X = Y(1); % Biomass, x
    P = Y(2); % Product, P
    S = Y(3); % Sustrate, S
    B = Y(4); % Cellobiose, B
    %Growth Kinetic equations
    
    %Product formation kineticss
%b(1) = k2
%b(2) = k2G
%b(3) = kG
%b(4) = um
%b(5) = m
%b(6) = k3
%b(7) = YXG
%b(8) = k1
%b(9) = e
%b(10) = k1G
%b(11) = k1B
%b(12) = k4
    %Material balance equations
    C=40-(0.9*S)-(0.947*B)-(0.9*P/0.511)-1.137*(X-1.04);
    dXdt = (b(4)*X*S)/(b(3)+S);
    dPdt = (b(6)*b(4)*X*S)/(b(3)+S)*b(7);
    dSdt = ((b(1)*B)/(0.95*(1+(S/b(2)))))-((b(4)*X*S)/(b(3)+S)*b(7))-(b(5)*X);
    dBdt = (b(8)*C)/(0.947*(1+(S/b(10))+(B/b(11))))*(b(9)/(1+b(9)*(12)*t))-((b(1)*B)/(1+(S/b(2)))); 
    dYdt = [dXdt;dPdt;dSdt;dBdt];
end

[Alfred]

@Torsten

I have not tried using the initial data from the article. I will try that ASAP.

[Star Strider]

I am running your code in MATLAB Online (and have been for the past 30 minutes at this point).

I changed the ODE solver to ode23s, since most kinetic problems are ‘stiff’, with parameter amplitude variations over several orders-of-magnitude. That will not change the result, however it should considerably shorten the time that the code needs to run.

14 Jan 2025 at 22:50

[Alfred]

Okay. Good. So, is ther any progress? @Star Strider

[Star Strider]

Unfortuantely, no.

I’ve had it running all night (>12 hours) and there is no indication it’s converging.

I’m changing the solve call to:

[sol,optval] = solve(prob,initialGuess,MultiStart,MinNumStartPoints=100)

with the expectation that considering different parameeter sets might be more likely to find a solution.

This requires the Global Optimization Toolbox, however I have that so this will work for me. What you need just now are a working set of parameters.

.

[Alfred]

Okay. That's fine.

I am also suggesting using Lsqnonlin, because this is what is normally reported in literature in related to my work. But, I don't know how to use it. @Star Strider

[Star Strider]

This is working, after I tweaked it a bit. (I’m not used to using the problem-based approach, so I’m learning that as I proceed here. My usual approach is the one I used in the earlier example that I posted.)

I appended (horzcat) ‘Y0’ to ‘X0’ since in my experience, estimating the initial conditions as parameters works best. I changed the rest of the code to accommodate that.

So:

b = optimvar('b',16,"LowerBound",0,"UpperBound",Inf);

and:

initialGuess.b = [3.20 16.25 3.40 0.18 0.64 3.02 0.52 0.61 1.00 0.18 5.85 0.10 Y0];

and in ‘toODE’:

Y0 = b(13:16);

The current solve call:

opts = optimoptions('lsqnonlin', Display='iter');

[sol,optval] = solve(prob, initialGuess, MultiStart, Options=opts, MinNumStartPoints=100)

That seems to be working better.

.

[Alfred]

Okay. Great.

Please could you send all the codes here. I am getting errors when I tried adding it to mine. @Star Strider

[Star Strider]

This is what I am currently running:

%% Fitting ODE Parameters (Batch Fermentation)

% Load experimental data

file = websave('exodefit', 'https://www.mathworks.com/matlabcentral/answers/uploaded_files/1823244/exodefit.xlsx')

%Experimantal Biomass data

% Xexp = readmatrix('exodefit.xlsx','Sheet','2','Range','d5:d18');

% Pexp = readmatrix('exodefit.xlsx','Sheet','2','Range','e5:e18');

% Sexp = readmatrix('exodefit.xlsx','Sheet','2','Range','f5:f18');

% Bexp = readmatrix('exodefit.xlsx','Sheet','2','Range','g5:g18');

Xexp = readmatrix(file,'Sheet','2','Range','d5:d18');

Pexp = readmatrix(file,'Sheet','2','Range','e5:e18');

Sexp = readmatrix(file,'Sheet','2','Range','f5:f18');

Bexp = readmatrix(file,'Sheet','2','Range','g5:g18');

% Group all data into yariable Y

Yexp =[Xexp'; Pexp'; Sexp'; Bexp'];

%Exprimantal time data

% timeexp = readmatrix('exodefit.xlsx','Sheet','2','Range','c5:c18');

timeexp = readmatrix(file,'Sheet','2','Range','c5:c18');

% Set ODE solver

% Define initial condition for ODE solver

tic

%initial biomass concentration

Y0 = [1.04 0.00 0.78 0.00];

%%

% Define time range for ODE solver

%Fermentation time

tspan = linspace(0,96);

% Set optimization problem

% Declare b = matrix of paramters to be determined

% b= [k2 k2G kG um m k3 YXG k1 e k1G k1B k4]

% b = optimvar('b',12,"LowerBound",0,"UpperBound",100);

b = optimvar('b',16,"LowerBound",0,"UpperBound",Inf);

%%

% Convert function for ODE solving to optimization expression

%myfcn = solution from ODE solving

myfcn = fcn2optimexpr(@toODE,b,timeexp)

% myfcn = fcn2optimexpr(@toODE,b,timeexp,Y0)

%%

% Set sum of squared error as an objective function to be minimized

SSE = sum(sum((myfcn-Yexp).^2));

%%

% Set optimization problem

prob = optimproblem ("Description","Fit ODE parameters",'ObjectiveSense','min')

%%

% Objective function (to be minimized)

prob.Objective = SSE

%%

% Show structure of problem

show(prob)

%%

% Set initial guesses of parameters

%b = [k2 k2G kG um m k3 YXG]

initialGuess.b = [3.20 16.25 3.40 0.18 0.64 3.02 0.52 0.61 1.00 0.18 5.85 0.10 Y0];

% initialGuess.b = [3.20 16.25 3.40 0.18 0.64 3.02 0.52 0.61 1.00 0.18 5.85 0.10];

%b(1) = k2

%b(2) = k2G

%b(3) = kG

%b(4) = um

%b(5) = m

%b(6) = k3

%b(7) = YXG

%b(8) = k1

%b(9) = e

%b(10) = k1G

%b(11) = k1B

%b(12) = k4

%%

% Solve optimization problem

% [sol,optval] = solve(prob,initialGuess)

opts = optimoptions('lsqnonlin', Display='iter');

[sol,optval] = solve(prob, initialGuess, MultiStart, Options=opts, MinNumStartPoints=100)

%%

toc

% Extract the results

%Fitted Parameters

bfinal =sol.b

%Sum of square Error

SSEfinal = optval

% Plot the simulated data and experimental data

%Call ODE to solve an equation using Final Fitted Parameters (bfinal)

solYsim = ode23s(@(t,Y)fermenterODEs(t,Y,bfinal),tspan)

%Evaluate the simulated results at specified tspan

Ysim = deval(solYsim,tspan)

%%

% Plot graphs

figure

tiledlayout(4,1)

nexttile

plot(timeexp,Xexp,'*r',tspan,Ysim(1,:),'-b')

xlabel('Time (h)')

ylabel('X (g/L)')

grid

nexttile

plot(timeexp,Pexp','*r',tspan,Ysim(2,:),'-b')

xlabel('Time (h)')

ylabel('P (g/L)')

grid

nexttile

plot(timeexp,Sexp','*r',tspan,Ysim(3,:),'-b')

xlabel('Time (h)')

ylabel('S (g/L)')

grid

nexttile

plot(timeexp,Bexp','*r',tspan,Ysim(4,:),'-b')

xlabel('Time (h)')

ylabel('B (g/L)')

grid

% Set functions for ODE solver

% This function is for solving ODE

function solferment = toODE(b,tspan)

% function solferment = toODE(b,tspan,Y0)

%call ODE45

Y0 = b(13:16);

sol = ode23s(@(t,Y)fermenterODEs(t,Y,b),tspan,Y0);

%Evaluate solution using tspan

solferment = deval(sol,tspan);

end

%%

% This function is for writing equation

%b = [k2 k2G kG um m k3 YXG]

%b(1) = k2

%b(2) = k2G

%b(3) = kG

%b(4) = um

%b(5) = m

%b(6) = k3

%b(7) = YXG

%b(8) = k1

%b(9) = e

%b(10) = k1G

%b(11) = k1B

%b(12) = k4

function dYdt = fermenterODEs(t,Y,b)

X = Y(1); % Biomass, x

P = Y(2); % Product, P

S = Y(3); % Sustrate, S

B = Y(4); % Cellobiose, B

%Growth Kinetic equations

%Product formation kineticss

%b(1) = k2

%b(2) = k2G

%b(3) = kG

%b(4) = um

%b(5) = m

%b(6) = k3

%b(7) = YXG

%b(8) = k1

%b(9) = e

%b(10) = k1G

%b(11) = k1B

%b(12) = k4

%Material balance equations

C=40-(0.9*S)-(0.947*B)-(0.9*P/0.511)-1.137*(X-1.04);

dXdt = (b(4)*X*S)/(b(3)+S);

dPdt = (b(6)*b(4)*X*S)/(b(3)+S)*b(7);

dSdt = ((b(1)*B)/(0.95*(1+(S/b(2)))))-((b(4)*X*S)/(b(3)+S)*b(7))-(b(5)*X);

dBdt = (b(8)*C)/(0.947*(1+(S/b(10))+(B/b(11))))*(b(9)/(1+b(9)*(12)*t))-((b(1)*B)/(1+(S/b(2))));

dYdt = [dXdt;dPdt;dSdt;dBdt];

end

It runs without error, however it is still not converging (after >2 hours).

Using websave avoids having to download the data file.

NOTE — You need to have the Global Optimization Toolbox for this version.

.

[Alfred]

Thank you @Star Strider

I have the global optimization toolbox.

[Star Strider]

My pleasure!

O.K., good!

In that context, here is a version of the code I previously popsted that uses the genetic algorithm —

function C=kinetics(theta,t)

c0 = theta(7:10);

% c0=[1;0;0;0];

[T,Cv]=ode45(@DifEq,t,c0);

%

function dC=DifEq(t,c)

dcdt=zeros(4,1);

dcdt(1)=-theta(1).*c(1)-theta(2).*c(1);

dcdt(2)= theta(1).*c(1)+theta(4).*c(3)-theta(3).*c(2)-theta(5).*c(2);

dcdt(3)= theta(2).*c(1)+theta(3).*c(2)-theta(4).*c(3)+theta(6).*c(4);

dcdt(4)= theta(5).*c(2)-theta(6).*c(4);

dC=dcdt;

end

C=Cv;

end

t=[0.1

0.2

0.4

0.6

0.8

1

1.5

2

3

4

5

6];

c=[0.902 0.06997 0.02463 0.00218

0.8072 0.1353 0.0482 0.008192

0.6757 0.2123 0.0864 0.0289

0.5569 0.2789 0.1063 0.06233

0.4297 0.3292 0.1476 0.09756

0.3774 0.3457 0.1485 0.1255

0.2149 0.3486 0.1821 0.2526

0.141 0.3254 0.194 0.3401

0.04921 0.2445 0.1742 0.5277

0.0178 0.1728 0.1732 0.6323

0.006431 0.1091 0.1137 0.7702

0.002595 0.08301 0.08224 0.835];

% theta0=[1;1;1;1;1;1];

% % [theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);

ftns = @(theta) norm(c-kinetics(theta,t));

PopSz = 500;

Parms = 10;

optsAns = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10); % Options Structure For 'Answers' Problems

opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10, 'PlotFcn',@gaplotbestf, 'PlotInterval',1); % Usual Options Call

optshf = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10, 'HybridFcn',@fmincon, 'PlotFcn',@gaplotbestf, 'PlotInterval',1); % With 'HybridFcn'

t0 = clock;

fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)

[theta,fval,exitflag,output,population,scores] = 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('Fitness value at convergence = %.4f\nGenerations \t\t\t\t = %d\n\n',fval,output.generations)

fprintf(1,'\tRate Constants:\n')

for k1 = 1:length(theta)

fprintf(1, '\t\tTheta(%2d) = %8.5f\n', k1, theta(k1))

end

tv = linspace(min(t), max(t));

Cfit = kinetics(theta, tv);

figure

hd = plot(t, c, 'p');

for k1 = 1:size(c,2)

CV(k1,:) = hd(k1).Color;

hd(k1).MarkerFaceColor = CV(k1,:);

end

hold on

hlp = plot(tv, Cfit);

for k1 = 1:size(c,2)

hlp(k1).Color = CV(k1,:);

end

hold off

grid

xlabel('Time')

ylabel('Concentration')

legend(hlp, compose('C_%d',1:size(c,2)), 'Location','N')

Consider this option (or soemetthing likee it), since I still can’t get the current code to converge after several hours.

.

[Alfred]

ok.

I will try it and get back to you.

Thanks @Star Strider

[Star Strider]

Thank you!

Something seems to be wrong, since it should have converged by now. The code doesn’t seem to be calling lsqnonlin. I have it set to report the iteration results, and none have appeared since the first set, just after it started. (I am not certain how the problem-oriented optimisation works, so this may be normal. The approach I am used to using should be calling the optimisation function repeatedly.)

[Torsten]

@Alfred

I really wonder why you don't apply your code first on a problem that has been tested and published. For me, this would be the reasonable way to proceed because as long as the code doesn't work for the published problem, it will not work for your problem, either.

[Alfred]

Yes @Torsten, I thought likewise but I don't know how to go about it as I don't have their experimental data.

You can be of help on how I should go about it.

Thanks.

[Torsten]

Why do you need experimental data ? I just meant to run your ODE model with the parameters and initial conditions used in the publication and compare the resulting curves for X, G, B and E. This would at least indicate that your physical model is correctly implemented.

[Alfred]

Okay, I will try that.

Thanks @Torsten

[Alfred]

https://www.mathworks.com/matlabcentral/answers/2145384-how-to-fit-kinetic-model-to-estimate-kinetic-parameter-from-experimental-data#comment_3322164

@Torsten and @Star Strider

I am using the same model in the publication.

I really don't know what to do. Can you guide me. I watched many videos on Youtube, but I am not getting it.

I have a few weeks to complete my work and this is the only challenge.

[Alfred]

If there is any way aside the code I write, kindly help. All models are from that paper.

[Torsten]

Please insert your code together with the initial conditions that reproduces the curves from the paper.

I don't know if it will make your subsequent fitting successful, but an error-free model is the prerequisite.

[Alfred]

@Star Strider and @Star Strider

would like to express my gratitude to everyone for your assistance on this journey.

I have just one week to resolve this issue, and any suitable approach to achieve this would be greatly appreciated.

I recently came across a work similar to mine, with the same objective: performing parameter estimation and curve fitting. You can find the link here: https://www.mathworks.com/matlabcentral/answers/314965-parameter-estimation-for-a-system-of-differential-equations.

I have edited the code from the link (Igor_Moura.m) above to accommodate my variables, but it has been running for over 5 hours without completing.

Your feedback and suggestions on this matter would be immensely helpful.

Please find the attached .m file for reference. Thank you!

Igor_Moura()

Local minimum possible.
lsqcurvefit stopped because the size of the current step is less than
the value of the step size tolerance.

	Rate Constants:
		Theta(1) =  1.00000
		Theta(2) =  1.00000
		Theta(3) =  1.00000
		Theta(4) =  1.00000
		Theta(5) =  1.00000
		Theta(6) =  1.00000
		Theta(7) =  1.00000
		Theta(8) =  1.00000
		Theta(9) =  1.00000
		Theta(10) =  1.00000
		Theta(11) =  1.00000
		Theta(12) =  1.00000

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

function C=kinetics(theta,t)
c0=[0.78;0;0;1.04];
[T,Cv]=ode15s(@DifEq,t,c0);
%
    function dC=DifEq(t,c)
    dcdt=zeros(4,1);
    C=40-(0.9*c(1))-(0.947*c(2))-(0.9*c(3)/0.511)-1.137*(c(4)-1.04);
    dcdt(1)=(theta(4)*c(4)*c(1))/(theta(3)+c(1));
    dcdt(2)= (theta(6)*theta(4)*c(4)*c(1))/((theta(3)+c(1))*theta(7));
    dcdt(3)= ((theta(1)*c(2))/(0.95*(1+(c(1)/theta(2)))))-((theta(4)*c(4)*c(1))/(theta(3)+c(1))*theta(7))-(theta(5)*c(4));
    dcdt(4)= (theta(8)*C)/(0.947*(1+(c(1)/theta(9))+(c(2)/theta(10))))*(4/(1+(4*theta(11)*t)))-((theta(1)*c(2))/(1+(c(1)/theta(2)))); 
    dC=dcdt;
    end
C=Cv;
end


t=[0
2
4
6
8
12
18
24
36
48
60
72
84
96];

c=[0.78	0.00	0.00	1.04
    7.44	1.43	3.19	1.05
11.75	2.88	3.46	1.31
15.99	3.09	3.56	1.32
20.03	1.32	4.68	1.28
20.75	2.09	7.99	1.29
19.62	1.51	10.36	1.27
17.87	0.25	9.89	2.29
13.75	1.09	12.24	2.75
6.26	0.28	12.98	2.94
5.17	0.46	13.74	2.88
5.35	0.23	17.86	2.12
5.07	1.01	19.22	1.04
5.07	0.71	17.99	1.01]; 
theta0=[1;1;1;1;1;1;1;1;1;1;1;1];

[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);

fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end

tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);

figure(1)
plot(t, c, 'p')
hold on
hlp = plot(tv, Cfit);
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'C_1(t)', 'C_2(t)', 'C_3(t)', 'C_4(t)', 'Location','N')

end

[Star Strider]

@Alfred — That is the essence of the code I posted here, both as an example and with your differential equation system. (I wrote that Answer.)

A version of it using the genetic algorithm (ga) function —

function C=kinetics(theta,t)
c0 = theta(7:10);
% c0=[1;0;0;0];
[T,Cv]=ode45(@DifEq,t,c0);
%
    function dC=DifEq(t,c)
    dcdt=zeros(4,1);
    dcdt(1)=-theta(1).*c(1)-theta(2).*c(1);
    dcdt(2)= theta(1).*c(1)+theta(4).*c(3)-theta(3).*c(2)-theta(5).*c(2);
    dcdt(3)= theta(2).*c(1)+theta(3).*c(2)-theta(4).*c(3)+theta(6).*c(4);
    dcdt(4)= theta(5).*c(2)-theta(6).*c(4);
    dC=dcdt;
    end
C=Cv;
end



t=[0.1
0.2
0.4
0.6
0.8
1
1.5
2
3
4
5
6];

c=[0.902	0.06997	0.02463	0.00218
0.8072	0.1353	0.0482	0.008192
0.6757	0.2123	0.0864	0.0289
0.5569	0.2789	0.1063	0.06233
0.4297	0.3292	0.1476	0.09756
0.3774	0.3457	0.1485	0.1255
0.2149	0.3486	0.1821	0.2526
0.141	0.3254	0.194	0.3401
0.04921	0.2445	0.1742	0.5277
0.0178	0.1728	0.1732	0.6323
0.006431	0.1091	0.1137	0.7702
0.002595	0.08301	0.08224	0.835]; 
% theta0=[1;1;1;1;1;1];

% % [theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
ftns = @(theta) norm(c-kinetics(theta,t));
PopSz = 250;
Parms = 10;

optsAns = optimoptions(@ga, 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10);                                                                  % Options Structure For 'Answers' Problems
opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);
optshf = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10, 'HybridFcn',@fmincon, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);   % With 'HybridFcn'

t0 = clock;
fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)

Start Time: 2025-01-27 20:15:51.4325

[theta,fval,exitflag,output,population,scores] = ga(ftns, Parms, [],[],[],[],zeros(Parms,1),Inf(Parms,1),[],[],optsAns);

ga stopped because the average change in the fitness value is less than options.FunctionTolerance.

t1 = clock;
fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)

Stop Time: 2025-01-27 20:17:00.3764

GA_Time = etime(t1,t0)

GA_Time = 68.9439

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)

Elapsed Time:   6.894389400000000E+01		00:01:08.943

fprintf('Fitness value at convergence = %.4f\nGenerations \t\t\t\t = %d\n\n',fval,output.generations)

Fitness value at convergence = 0.1923
Generations 				 = 1212

fprintf(1,'\tRate Constants:\n')

	Rate Constants:

for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%2d) = %8.5f\n', k1, theta(k1))
end

		Theta( 1) =  0.10701
		Theta( 2) =  0.90423
		Theta( 3) =  3.12931
		Theta( 4) =  7.44808
		Theta( 5) =  0.94748
		Theta( 6) =  0.23174
		Theta( 7) =  0.94321
		Theta( 8) =  0.06507
		Theta( 9) =  0.08379
		Theta(10) =  0.00000

tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);

figure
hd = plot(t, c, 'p');
for k1 = 1:size(c,2)
    CV(k1,:) = hd(k1).Color;
    hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Cfit);
for k1 = 1:size(c,2)
    hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, compose('C_%d',1:size(c,2)), 'Location','N')

Tightening the tolerances will give better results. You can use the ga results as initial estimates to a gradient descent function (such as lsqcurvefit).

.

[Alfred]

@Star Strider and @Torsten

The plot should like the file attached, that is when I just plot the raw data.

Look at the corrected codes...

@Star Strider, using the ga seems to be good but the ranges of time and temperature is not the same as the data.

function Igor_edit


function C=kinetics(theta,t)
c0=[0.78;0;0;1.04];
[T,Cv]=ode45(@DifEq,t,c0);
%
    function dC=DifEq(t,c)
    dcdt=zeros(4,1);
    C=40-(0.9*c(1))-(0.947*c(2))-(0.9*c(3)/0.511)-1.137*(c(4)-1.04);
    dcdt(1)= ((theta(1)*c(2))/(0.95*(1+(c(1)/theta(2)))))-((theta(4)*c(4)*c(1))/(theta(3)+c(1))*theta(7))-(theta(5)*c(4));
    dcdt(2)= (theta(8)*C)/(0.947*(1+(c(1)/theta(9))+(c(2)/theta(10))))*(4/(1+(4*theta(11)*t)))-((theta(1)*c(2))/(1+(c(1)/theta(2)))); 
    dcdt(3)= (theta(6)*theta(4)*c(4)*c(1))/((theta(3)+c(1))*theta(7));
    dcdt(4)= (theta(4)*c(4)*c(1))/(theta(3)+c(1)); 
    dC=dcdt;
    end
C=Cv;
end


t=[0
    2
4
6
8
12
18
24
36
48
60
72
84
96];

c=[0.78	0.00	0.00	1.04
    7.44	1.43	3.19	1.05
11.75	2.88	3.46	1.31
15.99	3.09	3.56	1.32
20.03	1.32	4.68	1.28
20.75	2.09	7.99	1.29
19.62	1.51	10.36	1.27
17.87	0.25	9.89	2.29
13.75	1.09	12.24	2.75
6.26	0.28	12.98	2.94
5.17	0.46	13.74	2.88
5.35	0.23	17.86	2.12
5.07	1.01	19.22	1.04
5.07	0.71	17.99	1.01]; 
theta0=[1;1;1;1;1;1;1;1;1;1;1;1];

[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);

fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end

tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);

figure(1)
plot(t, c, 'p')
hold on
hlp = plot(tv, Cfit);
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'C_1(t)', 'C_2(t)', 'C_3(t)', 'C_4(t)', 'Location','N')

end

[Star Strider]

I have been wotrking with this for the past four hours, and I cannot get it to converge. I even tried with the genetic algorithm —

% % % Alfred_2025_01_28.m

clear all
close all

% function Igor_edit


function C=kinetics(theta,t)
% c0=[0.78;0;0;1.04];
c0 = theta(12:15);
[T,Cv]=ode45(@DifEq,t,c0);
%
    function dC=DifEq(t,c)
    dcdt=zeros(4,1);
    C=40-(0.9*c(1))-(0.947*c(2))-(0.9*c(3)/0.511)-1.137*(c(4)-1.04);
    dcdt(1)= ((theta(1)*c(2))/(0.95*(1+(c(1)/theta(2)))))-((theta(4)*c(4)*c(1))/(theta(3)+c(1))*theta(7))-(theta(5)*c(4));
    dcdt(2)= (theta(8)*C)/(0.947*(1+(c(1)/theta(9))+(c(2)/theta(10))))*(4/(1+(4*theta(11)*t)))-((theta(1)*c(2))/(1+(c(1)/theta(2)))); 
    dcdt(3)= (theta(6)*theta(4)*c(4)*c(1))/((theta(3)+c(1))*theta(7));
    dcdt(4)= (theta(4)*c(4)*c(1))/(theta(3)+c(1)); 
    dC=dcdt;
    end
C=Cv;
end


t = [ 0
      2
      4
      6
      8
     12
     18
     24
     36
     48
     60
     72
     84
     96];

c=[0.78	0.00	0.00	1.04
   7.44	1.43	3.19	1.05    
   11.75	2.88	3.46	1.31
   15.99	3.09	3.56	.32
   20.03	1.32	4.68	1.28
   20.75	2.09	7.99	1.29
   19.62	1.51	10.36	1.27
   17.87	0.25	9.89	2.29
   13.75	1.09	12.24	2.75
   6.26	    0.28	12.98	2.94
   5.17	    0.46	13.74	2.88
   5.35	    0.23	17.86	2.12
   5.07	    1.01	19.22	1.04
   5.07	    0.71	17.99	1.01]; 
% theta0=[1;1;1;1;1;1;1;1;1;1;1;1];
% % theta0 = [rand(11,1); c(1,:).']
% % 
% % [theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);

ftns = @(theta) norm(c-kinetics(theta,t));
PopSz = 100;
Parms = 15;

optsAns = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10);                                                                  % Options Structure For 'Answers' Problems
opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-1, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);
optshf = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10, 'HybridFcn',@fmincon, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);   % With 'HybridFcn'

opts.Display='diagnose';

t0 = clock;
fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)
[theta,fval,exitflag,output,population,scores] = 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('Fitness value at convergence = %.4f\nGenerations \t\t\t\t = %d\n\n',fval,output.generations)

fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end

tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);

figure(1)
plot(t, c, 'p')
hold on
hlp = plot(tv, Cfit);
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'C_1(t)', 'C_2(t)', 'C_3(t)', 'C_4(t)', 'Location','N')

% end

This is not converging regardless of how I define the 'InitialPopulationMatrix' values. It is still on its first generation after two hours (third try). I even reduced the population size to 100 to speed it up, and that is not working either. The model is just not fitting the data.

I am not certainn what to suggest at this point, other than running your model with the published parameters to see if it reproduces the curves in the paper with them. When it does, then run it initially with your data with the genetic algorithm, and use those parameter estimiates as the starting values to the lsqcurvefit version of this code to refine the estimates.

.

[Alfred]

Okay. That will be okay to run with the published parameters. I don't know how to do that though. Your help would be gladly appreciated as always.

Thank you @Star Strider

Answer 2

[Torsten]

I replaced

C = b(11)-(0.9*y(2))-(0.947*y(3))-(0.9*y(4)/0.511)-(1.137*(y(1)-1.04));

by

C = b(1)-(0.9*y(2))-(0.947*y(3))-(0.9*y(4)/0.511)-(1.137*(y(1)-1.04));

in your code. b(11) does not exist.

And I was not sure what to insert for dx/dt in

rG = ((1/b(8))*(dx/dt))-(b(9)*y(1));

I used

rG = 1/b(8)*u-b(9)*y(1);

Maybe it must be

rG = 1/b(8)*u*X-b(9)*y(1);

The ODE integrator has problems with your differential equations at t = 6.19.

And to restrict all your parameters to be between 0 and 72 is quite strange - especially because tspan ends at 72.

%Fermentation data

Xdata = readmatrix('Batch1.xlsx','Sheet','sheet1','Range','e2:e13');
Gdata = readmatrix('Batch1.xlsx','Sheet','sheet1','Range','b2:b13');
Bdata = readmatrix('Batch1.xlsx','Sheet','sheet1','Range','c2:c13');
Edata = readmatrix('Batch1.xlsx','Sheet','sheet1','Range','d2:d13');


% Group all data into yariable y
yinv =[Xdata'; Gdata'; Bdata'; Edata'];


%Data for time
timeex = readmatrix('Batch1.xlsx','Sheet','sheet1','Range','a2:a13');

%Set ODE solver
% Initial cell biomass concentration (Initial condition)
y0=[1.04 0.78 0.00 0.00];

%Fermentation time
tspan = linspace(0,72);

%Set optimization problem
%Let b = matrix of paramters to be determined
% b= [um ks k1 K1G K1B k2 K2G YXG m k3]
%b = optimvar('b',10,"LowerBound",0,"UpperBound",72);



%Convert function for ODE solving to optimization expression 

%To use RtoODE in an objective function, convert the function to an
%optimization expression by using fcn2optimexpr. 
%myfcn = fcn2optimexpr(@toODE,b,timeex,y0);

%Express the objective function as the sum of squared differences between
%the ODE solution and the solution with true parameters.
%SSE = sum(sum((myfcn-yinv).^2));

%Create an optimization problem with the objective function SSE.
%prob = optimproblem ("Description","Fit ODE parameters",'ObjectiveSense','min');

%Objective function (to be minimized)
%prob.Objective = SSE;

%Show structure of problem
%show(prob)

%Solve Problem 
%To find the best-fitting parameters x, give an initial guess
%x0 for the solver and call solve. 
% Set initial guesses of parameters
b0 = [0.18 1.0 0.61 0.18 5.85 3.20 16.25 0.11 3.40 3.02];

%Solve optimization problem
%[sol,optval] = solve(prob,initialGuess);
[bfinal,SSEfinal] = lsqcurvefit(@(b,xdata)toODE(b,xdata,y0),b0,timeex,yinv)

Warning: Failure at t=6.192305e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.421085e-14) at time t.

Error using lsqcurvefit (line 312)
Function value and YDATA sizes are not equal.

%Extract the results
%Fitted Parameters
%bfinal =sol.b;

%Sum of square Error
%SSEfinal = optval;


%Plot the simulated data and experimental data
%Call ODE to solve an equation using Final Fitted Parameters (bfinal)
[~,ysim] = ode15s(@(t,y)batchferment(t,y,bfinal),tspan,y0);
%Evaluate the simulated results at specified tspan
%ysim = deval(solysim,tspan);

%Plot graphs
figure;
plot(timeex,Xdata,'*r',tspan,ysim(1,:),'-b')
legend('exp','sim')
xlabel('Time (h)')
ylabel('X (g/L)')

figure;
plot(timeex,Gdata','*r',tspan,ysim(2,:),'-b')
legend('exp','sim')
xlabel('Time (h)')
ylabel('G (g/L)')

figure;
plot(timeex,Bdata','*r',tspan,ysim(3,:),'-b')
legend('exp','sim')
xlabel('Time (h)')
ylabel('B (g/L)')

figure;
plot(timeex,Edata','*r',tspan,ysim(4,:),'-b')
legend('exp','sim')
xlabel('Time (h)')
ylabel('E (g/L)')

%Set functions for ODE solver for solving ODE
function solferment = toODE(b,timeex,y0)
    [~,solferment] = ode15s(@(t,y)batchferment(t,y,b),timeex,y0);
    %solferment = deval(sol,tspan);
end
%Equations for Batch 
%y(1) = X = Biomass Concentration (g/l)
%y(2) = G = Glucose Concentration (g/l)
%y(3) = B = Cellobiose Concentration (g/L)
%y(4) = E = Ethanol Concentration (g/l)

function dydt = batchferment(t,y,b)

   X=y(1) ; 
   G=y(2) ;
   B=y(3) ;   
   E=y(4) ;

    %%%Growth equations
    %dx/dt = (um*X*G)/(KG+G)
    u = (b(1)*y(1)*y(2))/(b(2)+y(2));
    %u = (b(1)*X*G)/(b(2)+G);
    %parameters
  um=b(1); % (1/h);
  ks=b(2); % (g/L);


    %%%Cellobiose equations
    C = b(1)-(0.9*y(2))-(0.947*y(3))-(0.9*y(4)/0.511)-(1.137*(y(1)-1.04));
     r1 = (b(3)*C)/(1+(y(2)/b(4))+(y(3)/b(5)));
     r2 = (b(6)*y(3))/(1+(y(2)/b(7)));
  k1=b(3) ;
  K1G=b(4) ;
  K1B=b(5) ;
  k2=b(6) ;
  K2G=b(7) ;

  %C = Cellulose concentration (g/l)
  %C = C0-(0.9*G)-(0.947*B)-(0.9*E/0.511)-(1.137*(X-X0))
  %r1 = (k1*C)/(1+(G/K1G)+(B/K1B))
  %r2 = (k2*B)/(1+(G/K2G)

 %C = C0-(0.9*G)-(0.947*B)-(0.9*E/0.511)-(1.137*(X-X0))
  %C = 83.65-(0.9*y(2))-(0.947*y(3))-(0.9*y(4)/0.511)-(1.137*(y(1)-1.04))
  
    
    %%% Glucose equations
    rG = 1/b(8)*u-b(9)*y(1);
  YXG=b(8) ;
  m=b(9) ;
  %rG = ((1/YXG)*(dX/dt))-(m*X)

    
    %%% Ethanol Concentration
      k3=b(10) ;
    Et = k3*(u*X)*(1/YXG);

  

    %Material balance equations
    dXdt = u*X;
    dGdt = (r2/0.95)-rG;
    dBdt = (r1/0.947)-r2;
    dEdt = Et;
    dydt = [dXdt;dGdt;dBdt;dEdt];
end

[Alfred]

Thanks @Torsten. I am still facing errors.

[Alfred]

@Torsten

I chose the tspan to end at 72 because the time in the works ends at 72 hours. Is there something I don't understand with tspan? Kindly explain to me.

Thanks