How to solve 'Index exceeds matrix dimensions.' error?

for t = 1:990
    
    %Kinetic model
    Qs(t) = Qsmax*Cs(t)/(Ks+Cs(t))*(1-exp(-t/td));      %the error start from here
    Qolim(t) = Qomax*Co(t)/(Ko+Co(t))*Ki/(Ki+Ce(t));
    Qslim(t) = ucr/Yoxxs;
    A = Qolim(t)/Yos;
    Qsox(t) = min(min(Qs(t),Qslim(t)),A);
    Qsred(t) = Qs(t)-Qsox(t);
    Qeup(t) = Qemax*(Ce(t)/(Ke+Ce(t)))*(Ki/(Ki+Ce(t)));
    B = (Qolim(t)-Qsox(t)*Yos)*Yeo;
    Qeox(t) = min(Qeup(t),B); 
    Qepr(t) = Qsred(t)*Yes;
    u(t) = (Qsox(t)*Yoxxs)+(Qsred(t)*Yredxs)+(Qeox(t)*Yxe);
    Qc(t) = (Qsox(t)*Yoxcs)+(Qsred(t)*Yredcs)+(Qeox(t)*Yce);
    Qo(t) = (Qsox(t)*Yos)+(Qeox(t)*Yeo);
    RQ(t) = Qc(t)/Qo(t); 
        
    F(t) = F(t)*exp(a-t);       
    
    %Dynamic model
    dCs(t) = (F(t)/60/V(t)*(So-Cs(t)))-(((u(t)/Yoxxs)+(Qepr(t)/Yes)+Qm)*Cx(t));
    dCo(t) = (-Qo(t)*Cx(t))+(kLao(t)/60*(Coo-Co(t)))-(F(t)/60/V(t)*Co(t));
    dCe(t) = ((Qepr-Qeox)*Cx(t))-(F(t)/60/V(t)*Ce(t));
    dCx(t) = (u(t)*Cx(t))-(F(t)/60/V(t)*Cx(t));
    dV(t) = F(t)/60;
    kLao(t) = (113*(Fa/60/AR)^0.25)/60;
    
    Cs(t) = Cs(t)+dCs*0.06;
    Co(t) = Co(t)+dCo*0.06;
    Ce(t) = Ce(t)+dCe*0.06;
    Cx(t) = Cx(t)+dCx*0.06;
    V(t) = V(t)+dV*0.06;
    
end

[WalterWhite]

Try using the . operator for the calculations maybe?

[Daniel Pollard]

Either Qs or Cs contains a number of elements where that number is less than 990. At some point in the iteration, it's trying to call Qs(t) or Cs(t) (or both) and can't.

We can't fix the error because you haven't posted enough code.

[Bobby Fischer]

There's a bunch of variables which aren't predefined. In order to solve your problem you should give these too, Christine.

[Christine King]

%Variable declaration (initial condition)

Cx = 15; %Biomass concentration

Cs = 7; %Glucose concentration

V = 50; %Volume

So = 325; %Feed concentration

tf = 1:0.06:990; %final time 16.5hr = 990min, save data every 3.6s

Vfer = 100; %Volume of fermentator

Co = 7.54; %Oxygen concentration

Ce = 0; %Ethanol concentration

F = 0; %Feed rate

Fa = 100; %Air feed rate

td = 1; %Time delay

kLao = 0; %total volumetric mass transfer coefficient

a = 0

%Parameters in fermentation model

Ke = 0.1;

Ko = 9.6e-5;

Ki = 3.5;

Ks = 0.612;

Yoxxs = 0.585;

Yredxs = 0.05;

Yos = 0.3857;

Yeo = 1.1236;

Yes = 0.4859;

Yxe = 0.7187;

Yoxcs = 0.5744;

Yredcs = 0.462;

Yce = 0.645;

Qemax = 3.967e-3;

Qomax = 4.25e-3;

Qsmax = 0.04905;

Qm = 5e-4;

ucr = 3.5e-3;

So = 325;

Coo = 0.006; %Co*

AR = 12.56;

%Initial condition

dV = 0; %dV/dt

dCs = 0; %dCs/dt

dCo = 0; %dCo/dt

dCe = 0; %dCe/dt

dCx = 0; %dCx/dt

%Fermentation process

for t = 0:length(tf)

%Kinetic model

Qs(t+1) = Qsmax*Cs(t+1)/(Ks+Cs(t+1))*(1-exp(-t/td));

Qolim(t+1) = Qomax*Co(t+1)/(Ko+Co(t+1))*Ki/(Ki+Ce(t+1));

Qslim(t+1) = ucr/Yoxxs;

A = Qolim(t+1)/Yos;

Qsox(t+1) = min(min(Qs(t+1),Qslim(t+1)),A);

Qsred(t+1) = Qs(t+1)-Qsox(t+1);

Qeup(t+1) = Qemax*(Ce(t+1)/(Ke+Ce(t+1)))*(Ki/(Ki+Ce(t+1)));

B = (Qolim(t+1)-Qsox(t+1)*Yos)*Yeo;

Qeox(t+1) = min(Qeup(t+1),B);

Qepr(t+1) = Qsred(t+1)*Yes;

u(t+1) = (Qsox(t+1)*Yoxxs)+(Qsred(t+1)*Yredxs)+(Qeox(t+1)*Yxe);

Qc(t+1) = (Qsox(t+1)*Yoxcs)+(Qsred(t+1)*Yredcs)+(Qeox(t+1)*Yce);

Qo(t+1) = (Qsox(t+1)*Yos)+(Qeox(t+1)*Yeo);

RQ(t+1) = Qc(t+1)/Qo(t+1);

F(t+1) = F(t+1)*exp(a-t);

%Dynamic model

dCs(t+1) = (F(t+1)/60/V(t+1)*(So-Cs(t+1)))-(((u(t+1)/Yoxxs)+(Qepr(t+1)/Yes)+Qm)*Cx(t+1));

dCo(t+1) = (-Qo(t+1)*Cx(t+1))+(kLao(t+1)/60*(Coo-Co(t+1)))-(F(t+1)/60/V(t+1)*Co(t+1));

dCe(t+1) = ((Qepr-Qeox)*Cx(t+1))-(F(t+1)/60/V(t+1)*Ce(t+1));

dCx(t+1) = (u(t+1)*Cx(t+1))-(F(t+1)/60/V(t+1)*Cx(t+1));

dV(t+1) = F(t+1)/60;

kLao(t+1) = (113*(Fa/60/AR)^0.25)/60;

Cs(t+1) = Cs(t)+dCs*0.06;

Co(t+1) = Co(t)+dCo*0.06;

Ce(t+1) = Ce(t)+dCe*0.06;

Cx(t+1) = Cx(t)+dCx*0.06;

V(t+1) = V(t)+dV*0.06;

end

Sorry, this is the full code I had written with all the parameters. Can help me solving this issue?

Answer 1

[Bobby Fischer]

Hi, there was a bunch of things to correct. Don't know if the result still has some value because I don't know what I was working with.

clear
%Variable declaration (initial condition)
Cx(1) = 15;      %Biomass concentration
Cs(1) = 7;       %Glucose concentration
V(1) = 50;       %Volume
So = 325;     %Feed concentration
tf = 1:0.06:990;     %final time 16.5hr = 990min, save data every 3.6s
Vfer = 100;   %Volume of fermentator
Co(1) = 7.54;       %Oxygen concentration
Ce(1) = 0;       %Ethanol concentration
F(1)= 0;        %Feed rate %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Fa = 100;       %Air feed rate
td = 1;       %Time delay
kLao(1) = 0;       %total volumetric mass transfer coefficient
a = 0;
%Parameters in fermentation model
Ke = 0.1;
Ko = 9.6e-5;
Ki = 3.5;
Ks = 0.612;
Yoxxs = 0.585;
Yredxs = 0.05;
Yos = 0.3857;
Yeo = 1.1236;
Yes = 0.4859;
Yxe = 0.7187;
Yoxcs = 0.5744;
Yredcs = 0.462;
Yce = 0.645;
Qemax = 3.967e-3;
Qomax = 4.25e-3;
Qsmax = 0.04905;
Qm = 5e-4;
ucr = 3.5e-3;
So = 325;
Coo = 0.006;    %Co*
AR = 12.56;
%Initial condition
dV = 0;     %dV/dt
dCs(1) = 0;    %dCs/dt
dCo(1) = 0;    %dCo/dt
dCe(1) = 0;    %dCe/dt
dCx(1) = 0;    %dCx/dt
%Fermentation process
for t = 1:length(tf)-1
    
    %Kinetic model
    Qs(t+1) = Qsmax*Cs(t)/(Ks+Cs(t))*(1-exp(-t/td));
    Qolim(t+1) = Qomax*Co(t)/(Ko+Co(t))*Ki/(Ki+Ce(t));
    Qslim(t+1) = ucr/Yoxxs;
    A = Qolim(t+1)/Yos;
    Qsox(t+1) = min(min(Qs(t+1),Qslim(t+1)),A);
    Qsred(t+1) = Qs(t+1)-Qsox(t+1);
    Qeup(t+1) = Qemax*(Ce(t)/(Ke+Ce(t)))*(Ki/(Ki+Ce(t)));
    B = (Qolim(t+1)-Qsox(t+1)*Yos)*Yeo;
    Qeox(t+1) = min(Qeup(t+1),B);
    Qepr(t+1) = Qsred(t+1)*Yes;
    u(t+1) = (Qsox(t+1)*Yoxxs)+(Qsred(t+1)*Yredxs)+(Qeox(t+1)*Yxe);
    Qc(t+1) = (Qsox(t+1)*Yoxcs)+(Qsred(t+1)*Yredcs)+(Qeox(t+1)*Yce);
    Qo(t+1) = (Qsox(t+1)*Yos)+(Qeox(t+1)*Yeo);
    RQ(t+1) = Qc(t+1)/Qo(t+1);
    
    F(t+1) = F(t)*exp(a-t); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %Dynamic model
    dCs(t+1) = (F(t+1)/60/V(t)*(So-Cs(t)))-(((u(t+1)/Yoxxs)+(Qepr(t+1)/Yes)+Qm)*Cx(t));
    dCo(t+1) = (-Qo(t+1)*Cx(t))+(kLao(t)/60*(Coo-Co(t)))-(F(t+1)/60/V(t)*Co(t));
    dCe(t+1) = ((Qepr(t+1)-Qeox(t+1))*Cx(t))-(F(t+1)/60/V(t)*Ce(t)); %%%%%%%%%%%%%%%%%%%%%%%
    dCx(t+1) = (u(t+1)*Cx(t))-(F(t+1)/60/V(t)*Cx(t));
    dV(t+1) = F(t+1)/60;
    kLao(t+1) = (113*(Fa/60/AR)^0.25)/60;
    
    Cs(t+1) = Cs(t)+dCs(t+1)*0.06;
    Co(t+1) = Co(t)+dCo(t+1)*0.06;
    Ce(t+1) = Ce(t)+dCe(t+1)*0.06;
    Cx(t+1) = Cx(t)+dCx(t+1)*0.06;
    V(t+1) = V(t)+dV(t+1)*0.06;
    
end

whos
figure(1)
clf
hold on
plot(tf,Ce,'b')
plot(tf,Co,'r')
plot(tf,Cs,'g')
plot(tf,Cx,'y')
legend('Ce','Co','Cs','Cx')
figure(2)
clf
hold on
plot(tf,Qs,'b')
plot(tf,Qslim,'y')
plot(tf,Qsox,'g')
plot(tf,Qsred,'r')

legend('Qs','Qslim','Qsox','Qsred')