How to fit kinetic model to estimate kinetic parameter from experimental data
Afficher commentaires plus anciens
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'.
>>
Réponse acceptée
Star Strider
le 15 Août 2024
Modifié(e) : Star Strider
le 15 Août 2024
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.
.
65 commentaires
Star Strider
le 15 Août 2024
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
le 26 Août 2024
Star Strider
le 26 Août 2024
Modifié(e) : Star Strider
le 26 Août 2024
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)
[theta,fval,exitflag,output,population,scores] = ga(ftns, Parms, [],[],[],[],zeros(Parms,1),Inf(Parms,1),[],[],optsAns);
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 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
le 26 Août 2024
Modifié(e) : Walter Roberson
le 14 Jan 2025
Alfred
le 17 Oct 2024
Déplacé(e) : Star Strider
le 17 Oct 2024
Alfred
le 17 Oct 2024
Déplacé(e) : Star Strider
le 17 Oct 2024
Star Strider
le 17 Oct 2024
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
le 3 Nov 2024
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')
Alfred
le 3 Nov 2024
Torsten
le 3 Nov 2024
It would help if you included the .m-file instead of the .mlx-file.
Star Strider
le 3 Nov 2024
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.
.
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
le 7 Nov 2024
Star Strider
le 7 Nov 2024
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
le 7 Nov 2024
Star Strider
le 7 Nov 2024
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.)
.
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
le 8 Nov 2024
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)
[theta,fval,exitflag,output,population,scores] = ga(ftns, Parms, [],[],[],[],zeros(Parms,1),Inf(Parms,1),[],[],optsAns);
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
.
Star Strider
le 26 Nov 2024
Having the paper helps.
I wiill look through it (that will take a while) and make appropriate suggestions.
Alfred
le 26 Nov 2024
Star Strider
le 3 Déc 2024
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.
, , , 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
le 11 Jan 2025
Star Strider
le 11 Jan 2025
As always, my (our) pleasure!
Alfred
le 13 Jan 2025
Torsten
le 13 Jan 2025
Could you include the .m file, not the .mlx file ?
Alfred
le 14 Jan 2025
Alfred
le 14 Jan 2025
Star Strider
le 14 Jan 2025
Modifié(e) : Star Strider
le 14 Jan 2025
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)
%%
% 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
le 14 Jan 2025
Star Strider
le 14 Jan 2025
As always, my pleasure!
I have to be away for a while. I will get back to it when II return.
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
le 14 Jan 2025
Alfred
le 14 Jan 2025
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
le 14 Jan 2025
Star Strider
le 14 Jan 2025
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
le 15 Jan 2025
Star Strider
le 15 Jan 2025
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
le 15 Jan 2025
Star Strider
le 15 Jan 2025
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
le 15 Jan 2025
Star Strider
le 15 Jan 2025
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
le 15 Jan 2025
Star Strider
le 15 Jan 2025
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
le 15 Jan 2025
Star Strider
le 15 Jan 2025
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.)
Alfred
le 15 Jan 2025
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
le 15 Jan 2025
Alfred
le 19 Jan 2025
Alfred
le 19 Jan 2025
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()
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
@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.)
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)
[theta,fval,exitflag,output,population,scores] = ga(ftns, Parms, [],[],[],[],zeros(Parms,1),Inf(Parms,1),[],[],optsAns);
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')
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
le 28 Jan 2025
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
le 28 Jan 2025
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
le 28 Jan 2025
Plus de réponses (1)
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)
%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
Catégories
En savoir plus sur Genetic Algorithm dans Centre d'aide et File Exchange
Produits
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
