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.
.
