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.

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