Clear Filters
Clear Filters

I am struggling to set a maximum value for a output for a function. I have nested my function in another but it is giving an error.

3 views (last 30 days)
I have the below script and main function. After nesting so as to introduce a limit for maxSnh4, it is giving an error. Any hint on how to restructure it would be highly appreciated.
%Script:
tSpan=[0 300];
y0=[1;1;1;1;1;1;1;1;1;1;1];
M1=eye(55);
options = odeset('RelTol',1e-20,'AbsTol',1e-10,'Mass',M1,'MassSingular', 'maybe','MStateDependence','strong','InitialStep',1e-10, 'MStateDependence','strong');
[tSol, ySol]=ode15s(@(t,y) SBR(t,y), tSpan, y0);
%plotting
figure (1)
plot(tSol,ySol(:,7),'-r'); hold on %Sse
plot(tSol,ySol(:,8),'-g'); hold on %Sno3
plot(tSol,ySol(:,9),'--r'); hold on %Sno2
plot(tSol,ySol(:,10),'-k'); hold on %nh4
plot(tSol,ySol(:,11),'-b'); hold on %So2
figure (2)
plot(tSol,ySol(:,1),'-r'); hold on %Xaob
plot(tSol,ySol(:,2),'-g'); hold on %Xnb
plot(tSol,ySol(:,3),'--r'); hold on %Xnsp
plot(tSol,ySol(:,4),'-k'); hold on %amx
plot(tSol,ySol(:,5),'-b'); hold on %Xhan
plot(tSol,ySol(:,6),'--g'); hold off %Xs
xlabel('Time (days)'); ylabel('Concentration (g/m^3)');
%Main function
% SBR simulation in MATLAB
function dy=SBR(t,y)
maxSnh4=200;
function fval=SBR1(t,y)
% Initialize arrays
Xaob=y(1);
Xnb = y(2);
Xnsp= y(3);
Xamx= y(4);
Xhan= y(5);
Xs= y(6);
Sse= y(7);
Sno3= y(8);
Sno2= y(9);
Snh4= y(10);
So2= y(11);
T=303;
%aob
O1=68/(0.00831*293*(273+T)); u1=0.054*exp(O1*(T-293)); Ko2aob=0.6; Knh4aob=2.4; Yaob=0.15; Kno3aob=0.5;
%nb
O2=44/(0.00831*293*(273+T)); u2=0.047*exp(O2*(T-293)); Ko2nb=2.2; Kno2nb=0.5*1.375; Ynb=0.041;Kno3nb=0.5;
%nsp
O3=44/(0.00831*293*(273+T)); u3=0.69*exp(O3*(T-293)); Ko2nsp=0.33; Kno2nsp=0.52; Ynsp=0.14;Kno3nsp=0.5;
%amx
O5=71/(0.00831*293*(273+T)); u5=0.0022*exp(O5*(T-293)); Kno2amx=0.05; Knh4amx=0.07; Yamx=0.159; Ko2amx=0.01; Kno3amx=0.5;
%han
O6=0.069;u6=0.5*7.2*exp((T-293)*0.069); KSsehan=2; Kno2han=0.5;Kno3han=0.5; YNo2han=0.3; Ko2han=0.2;
YNo3han=0.43; Yhaer=0.54;
%other constants
baob=0.0054*exp(O6*(T-293)); bnb=0.0025*exp(O6*(T-293));bnsp=bnb; bamx=0.00013*exp(O6*(T-293)); bhan=0.008*exp(O6*(T-293));
INBM=0.07; fI=0.1; namx=0.5; nhan=0.6;naob=0.5; nnb=0.5; nnsp=0.5;
KX=1; KH=0.125*exp(0.069*(T-293)); INXI=0.02;
% input parameters
Sseo=300; Sno3o=0.2; Qo=116.67; Xso=250;
Snh4o=650; Sno2o=0.2; So2o=0;
%oxygen transfer parameters
KaL=1920;
V=1400; SRT=100;
Xaobo=2; Xnbo=2; Xnspo=2; Xhano=2; Xamxo=2;
Vo=V-0.3*V; % volume in the reactor at the start of refilling
% SBR operation cycle
refilling_duration = 1.59; % Duration of refilling phase (min)
aeration_duration = 15; % Duration of aeration phase (min)
mixing_duration = 45; % Duration of mixing phase (min)
settling_duration = 30; % Duration of settling phase (min)
withdrawal_duration = 0.58; % Duration of withdrawal phase (min)
cycle_duration=refilling_duration+aeration_duration+mixing_duration+settling_duration+withdrawal_duration;
t0=0+cycle_duration;
t1=t0+15;
t2=t1+15;
t3=t2+45;
t4=t3+60;
t5=t4+15;
%refilling
while t>t0 && t<=t1
% input parameters
Sseo=300; Sno3o=0.2; Qo=220; Xso=250;
Snh4o=650; Sno2o=0.2; So2o=0; KaL=0;
Xaobo=2; Xnbo=2; Xnspo=2; Xhano=2; Xamxo=2;
break
end
%aeration
while t>t1 && t<=t2
% input parameters
Sseo=0; Sno3o=0; Qo=0; Xso=0; Snh4o=0; Sno2o=0; So2o=0; Xaobo=0; Xnbo=0; Xnspo=0; Xhano=0; Xamxo=0;
if Snh4<10 || Sno2>10 || So2>0.8
KaL=0;
elseif Snh4>10 && So2<=0.5
KaL=1920;
elseif Sno2<10 && So2<0.5
KaL=1920;
end
break
end
%mixing & settling
while t>t2 && t<=t4
% input parameters
Sseo=0; Sno3o=0; Qo=0; Xso=0; Snh4o=0; Sno2o=0.2; So2o=0; KaL=0; Xaobo=0; Xnbo=0; Xnspo=0; Xhano=0; Xamxo=0;
break
end
%withdrawal
while t>t4 && t<=t5
% input parameters
Sseo=0; Sno3o=0; Qo=600; Xso=0; Snh4o=0; Sno2o=0.2; So2o=0; KaL=0; Xaobo=0; Xnbo=0; Xnspo=0; Xhano=0; Xamxo=0;
break
end
fval(1)=((Vo*Xaobo-Vo*Xaob/SRT)/V + u1*(So2/(Ko2aob+So2))*(Snh4/(Knh4aob+Snh4))*Xaob - baob*(So2/(Ko2aob+So2))*Xaob - baob*naob*((Ko2aob/(Ko2aob+So2))*(Sno2+Sno3)/(Kno3aob+Sno2+Sno3))*Xaob);
fval(2)=((Vo*Xnbo-Vo*Xnb/SRT)/V + u2*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xnb - bnb*(So2/(Ko2nb+So2))*Xnb - bnb*nnb*((Ko2nb/(Ko2nb+So2))*(Sno2+Sno3)/(Kno3nb+Sno2+Sno3))*Xnb);
fval(3)=((Vo*Xnspo-Vo*Xnsp/SRT)/V + u3*(Sno2/(Kno2nsp+Sno2))*(So2/(Ko2nsp+So2))*Xnsp - bnsp*(So2/(Ko2nsp+So2))*Xnsp - bnsp*nnsp*((Ko2nsp/(Ko2nsp+So2))*(Sno2+Sno3)/(Kno3nsp+Sno2+Sno3))*Xnsp);
fval(4)=((Vo*Xamxo-Vo*Xamx/SRT)/V + u5*((Snh4/(Knh4amx+Snh4))*(Sno2/(Kno2amx+Sno2))*(Ko2amx/(Ko2amx+So2)))*Xamx-bamx*(So2/(Ko2amx+So2))*Xamx - bamx*namx*((Ko2amx/(Ko2amx+So2))*(Sno2+Sno3)/(Kno3amx+Sno2+Sno3))*Xamx);
fval(5)=((Vo*Xhano-Vo*Xhan/SRT)/V + (u6*(Sse/(KSsehan+Sse))*(Sno2/(Kno2han+Sno2))*(Ko2han/(Ko2han+So2)) + u6*(Sse/(KSsehan+Sse))*(Sno3/(Kno3han+Sno3))*(Ko2han/(Ko2han+So2))+u6*((Sse/(KSsehan+Sse))*(So2/(Ko2han+So2)))*Xhan)- bhan*(So2/(Ko2han+So2))*Xhan - bhan*nhan*((Ko2han/(Ko2han+So2))*(Sno2+Sno3)/(Kno3han+Sno2+Sno3))*Xhan);
fval(6)=((Qo*Xso-Qo*Xs)/V - KH*((Xs/Xhan)/(KX+(Xs/Xhan)))*Xhan);
fval(7)=((Qo*Sseo-Qo*Sse)/V - (u6*nhan*(1/YNo2han)*(Sse/(KSsehan+Sse))*(Sno2/(Kno2han+Sno2))*(Ko2han/(Ko2han+So2)))*Xhan -(u6*nhan*(1/YNo3han)*(Sse/(KSsehan+Sse))*(Sno3/(Kno3han+Sno3))*(Ko2han/(Ko2han+So2)))*Xhan - ((u6/Yhaer)*(Sse/(KSsehan+Sse))*(So2/(Ko2han+So2))*Xhan) + KH*((Xs/Xhan)/(KX+(Xs/Xhan)))*Xhan);
fval(8)=((Qo*Sno3o-Qo*Sno3)/V - (u6*nhan*((1-YNo3han)/(2.86*YNo3han))*(Sse/(KSsehan+Sse))*(Sno3/(Kno3han+Sno3))*(Ko2han/(Ko2han+So2)))*Xhan + (u5/1.14)*(Snh4/(Knh4amx+Snh4))*(Sno2/(Kno2amx+Sno2))*(Ko2amx/(Ko2amx+So2))*Xamx + (u2/Ynb*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xnb) + ((u3/Ynsp)*(Sno2/(Kno2nsp+Sno2))*(So2/(Ko2nsp+So2))*Xnsp) - ((1-fI)/2.86)*baob*naob*(Ko2aob/(Ko2aob+So2))*(Sno2+Sno3)/(Kno3aob+Sno2+Sno3)*Xaob - ((1-fI)/2.86)*bnb*nnb*(Ko2nb/(Ko2nb+So2))*(Sno2+Sno3)/(Kno3nb+Sno2+Sno3)*Xnb - ((1-fI)/2.86)*bnsp*nnsp*(Ko2nsp/(Ko2nsp+So2))*(Sno2+Sno3)/(Kno3nsp+Sno2+Sno3)*Xnsp - ((1-fI)/2.86)*bamx*namx*(Ko2amx/(Ko2amx+So2))*(Sno2+Sno3)/(Kno3amx+Sno2+Sno3)*Xamx-((1-fI)/2.86)*bhan*nhan*(Ko2han/(Ko2han+So2))*(Sno2+Sno3)/(Kno3han+Sno2+Sno3)*Xhan);
fval(9)=((Qo*Sno2o-Qo*Sno2)/V - (((1-YNo2han)/(1.71*YNo2han))*u6*nhan*(Sse/(KSsehan+Sse))*(Sno2/(Kno2han+Sno2))*(Ko2han/(Ko2han+So2))*Xhan) - (((u5/Yamx)+(u5/1.14))*(Snh4/(Knh4amx+Snh4))*(Sno2/(Kno2amx+Sno2))*(Ko2amx/(Ko2amx+So2)))*Xamx + ((u1/Yaob)*(So2/(Ko2aob+So2))*(Snh4/(Knh4aob+Snh4)))*Xaob - ((u2/Ynb)*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xnb) - ((u3/Ynsp)*(Sno2/(Kno2nsp+Sno2))*(So2/(Ko2nsp+So2)))*Xnsp);
fval(10)=((Qo*Snh4o-Qo*Snh4)/V - (u1*(INBM+1/Yaob)*(So2/(Ko2aob+So2))*(Snh4/(Knh4aob+Snh4)))*Xaob -(u5*(INBM+1/Yamx)*(Snh4/(Knh4amx+Snh4))*(Sno2/(Kno2amx+Sno2))*(Ko2amx/(Ko2amx+So2)))*Xamx - (u2*INBM*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2)))*Xnb -(u3*INBM*(Sno2/(Kno2nsp+Sno2))*(So2/(Ko2nsp+So2)))*Xnsp - INBM*u6*nhan*((Sse/(KSsehan+Sse))*(Sno2/(Kno2han+Sno2))*(Ko2han/(Ko2han+So2)))*Xhan - u6*INBM*nhan*((Sse/(KSsehan+Sse))*(Sno3/(Kno3han+Sno3))*(Ko2han/(Ko2han+So2)))*Xhan - u6*INBM*((Sse/(KSsehan+Sse))*(So2/(Ko2han+So2)))*Xhan + (INBM-(fI*INXI))*baob*naob*((Ko2aob/(Ko2aob+So2))*(Sno2+Sno3)/(Kno3aob+Sno2+Sno3))*Xaob + (INBM-(fI*INXI))*bnb*nnb*((Ko2nb/(Ko2nb+So2))*(Sno2+Sno3)/(Kno3nsp+Sno2+Sno3))*Xnb +(INBM-(fI*INXI))*bnsp*nnsp*((Ko2nsp/(Ko2nsp+So2))*(Sno2+Sno3)/(Kno3nsp+Sno2+Sno3))*Xnsp + (INBM-(fI*INXI))*bamx*namx*((Ko2amx/(Ko2amx+So2))*(Sno2+Sno3)/(Kno3amx+Sno2+Sno3))*Xamx + (INBM-(fI*INXI))*bhan*nhan*((Ko2han/(Ko2han+So2))*(Sno2+Sno3)/(Kno3han+Sno2+Sno3))*Xhan +(INBM-(fI*INXI))*baob*(So2/(Ko2aob+So2))*Xaob + (INBM-(fI*INXI))*bnb*(So2/(Ko2nb+So2))*Xnb +(INBM-(fI*INXI))*bnsp*(So2/(Ko2nsp+So2))*Xnsp + (INBM-(fI*INXI))*bamx*(So2/(Ko2amx+So2))*Xamx + (INBM-(fI*INXI))*bhan*(So2/(Ko2han+So2))*Xhan);
fval(11)=((Qo*So2o-Qo*So2)/V + KaL*(8-So2) -(((1-Yhaer)/Yhaer)*u6*(Sse/(KSsehan+Sse))*(So2/(Ko2han+So2)))*Xhan -(((3.43-Yaob)/Yaob)*u1*(So2/(Ko2aob+So2))*(Snh4/(Knh4aob+Snh4)))*Xaob - (((1.14-Ynb)/Ynb)*u2*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2)))*Xnb - (((1.14-Ynsp)/Ynsp)*u3*(Sno2/(Kno2nsp+Sno2))*(So2/(Ko2nsp+So2)))*Xnsp - (1-fI)*baob*(So2/(Ko2aob+So2))*Xaob - (1-fI)*bnb*(So2/(Ko2nb+So2))*Xnb -(1-fI)*bnsp*(So2/(Ko2nsp+So2))*Xnsp - (1-fI)*bamx*(So2/(Ko2amx+So2))*Xamx - (1-fI)*bhan*(So2/(Ko2han+So2))*Xhan);
if fval(10)>maxSnh4
dy=maxSnh4;
else
dy=fval(10);
end
end
end
  3 Comments
KIPROTICH KOSGEY
KIPROTICH KOSGEY on 13 Jul 2023
maxSnh4 is defined in row 3.
It is my first time coding nested functions and I am struggling to structure it properly hence the errors.
the error is: Output argument "dy" (and possibly others) not assigned a value in the execution with "SBR" function.
I just need to set fval(10) to <=200. Any hint would be greatly appreciated.
Stephen23
Stephen23 on 13 Jul 2023
"Output argument "dy" (and possibly others) not assigned a value in the execution with "SBR" function."
Because you never define the variable DY inside of the function SBR.
You also never call the nested function SBR1, so you might as well just delete it. Your code is exactly equivalent to this:
function dy=SBR(t,y)
maxSnh4=200;
end
which clearly does not define DY anywhere.

Sign in to comment.

Accepted Answer

Torsten
Torsten on 13 Jul 2023
Moved: Torsten on 13 Jul 2023
Instead of programming if-statements within the function routine and setting parameters depending on the actual time t, you should integrate up to t0, return control to the calling program and save results, reset parameters, integrate from t0 to t1, return control to the calling program and save results, reset parameters, integrate ... This is the usual procedure to avoid discontinuities in the ODE equations which most probably makes ode15s stop with an error message at time t ~ 107.
%Script:
tSpan=[0 300];
y0=[1;1;1;1;1;1;1;1;1;1;1];
M1=eye(55);
options = odeset('RelTol',1e-20,'AbsTol',1e-10,'Mass',M1,'MassSingular', 'maybe','MStateDependence','strong','InitialStep',1e-10, 'MStateDependence','strong');
[tSol, ySol]=ode15s(@(t,y) SBR(t,y), tSpan, y0);
Warning: Failure at t=1.071763e+02. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.273737e-13) at time t.
%plotting
figure (1)
plot(tSol,ySol(:,7),'-r'); hold on %Sse
plot(tSol,ySol(:,8),'-g'); hold on %Sno3
plot(tSol,ySol(:,9),'--r'); hold on %Sno2
plot(tSol,ySol(:,10),'-k'); hold on %nh4
plot(tSol,ySol(:,11),'-b'); hold on %So2
figure (2)
plot(tSol,ySol(:,1),'-r'); hold on %Xaob
plot(tSol,ySol(:,2),'-g'); hold on %Xnb
plot(tSol,ySol(:,3),'--r'); hold on %Xnsp
plot(tSol,ySol(:,4),'-k'); hold on %amx
plot(tSol,ySol(:,5),'-b'); hold on %Xhan
plot(tSol,ySol(:,6),'--g'); hold off %Xs
xlabel('Time (days)'); ylabel('Concentration (g/m^3)');
%Main function
% SBR simulation in MATLAB
function dy=SBR(t,y)
maxSnh4=200;
%function fval=SBR1(t,y)
% Initialize arrays
Xaob=y(1);
Xnb = y(2);
Xnsp= y(3);
Xamx= y(4);
Xhan= y(5);
Xs= y(6);
Sse= y(7);
Sno3= y(8);
Sno2= y(9);
Snh4= y(10);
So2= y(11);
T=303;
%aob
O1=68/(0.00831*293*(273+T)); u1=0.054*exp(O1*(T-293)); Ko2aob=0.6; Knh4aob=2.4; Yaob=0.15; Kno3aob=0.5;
%nb
O2=44/(0.00831*293*(273+T)); u2=0.047*exp(O2*(T-293)); Ko2nb=2.2; Kno2nb=0.5*1.375; Ynb=0.041;Kno3nb=0.5;
%nsp
O3=44/(0.00831*293*(273+T)); u3=0.69*exp(O3*(T-293)); Ko2nsp=0.33; Kno2nsp=0.52; Ynsp=0.14;Kno3nsp=0.5;
%amx
O5=71/(0.00831*293*(273+T)); u5=0.0022*exp(O5*(T-293)); Kno2amx=0.05; Knh4amx=0.07; Yamx=0.159; Ko2amx=0.01; Kno3amx=0.5;
%han
O6=0.069;u6=0.5*7.2*exp((T-293)*0.069); KSsehan=2; Kno2han=0.5;Kno3han=0.5; YNo2han=0.3; Ko2han=0.2;
YNo3han=0.43; Yhaer=0.54;
%other constants
baob=0.0054*exp(O6*(T-293)); bnb=0.0025*exp(O6*(T-293));bnsp=bnb; bamx=0.00013*exp(O6*(T-293)); bhan=0.008*exp(O6*(T-293));
INBM=0.07; fI=0.1; namx=0.5; nhan=0.6;naob=0.5; nnb=0.5; nnsp=0.5;
KX=1; KH=0.125*exp(0.069*(T-293)); INXI=0.02;
% input parameters
Sseo=300; Sno3o=0.2; Qo=116.67; Xso=250;
Snh4o=650; Sno2o=0.2; So2o=0;
%oxygen transfer parameters
KaL=1920;
V=1400; SRT=100;
Xaobo=2; Xnbo=2; Xnspo=2; Xhano=2; Xamxo=2;
Vo=V-0.3*V; % volume in the reactor at the start of refilling
% SBR operation cycle
refilling_duration = 1.59; % Duration of refilling phase (min)
aeration_duration = 15; % Duration of aeration phase (min)
mixing_duration = 45; % Duration of mixing phase (min)
settling_duration = 30; % Duration of settling phase (min)
withdrawal_duration = 0.58; % Duration of withdrawal phase (min)
cycle_duration=refilling_duration+aeration_duration+mixing_duration+settling_duration+withdrawal_duration;
t0=0+cycle_duration;
t1=t0+15;
t2=t1+15;
t3=t2+45;
t4=t3+60;
t5=t4+15;
%refilling
%while t>t0 && t<=t1
if t>t0 && t<=t1
% input parameters
Sseo=300; Sno3o=0.2; Qo=220; Xso=250;
Snh4o=650; Sno2o=0.2; So2o=0; KaL=0;
Xaobo=2; Xnbo=2; Xnspo=2; Xhano=2; Xamxo=2;
%break
end
%aeration
%while t>t1 && t<=t2
if t>t1 && t<=t2
% input parameters
Sseo=0; Sno3o=0; Qo=0; Xso=0; Snh4o=0; Sno2o=0; So2o=0; Xaobo=0; Xnbo=0; Xnspo=0; Xhano=0; Xamxo=0;
if Snh4<10 || Sno2>10 || So2>0.8
KaL=0;
elseif Snh4>10 && So2<=0.5
KaL=1920;
elseif Sno2<10 && So2<0.5
KaL=1920;
end
%break
end
%mixing & settling
%while t>t2 && t<=t4
if t>t2 && t<=t4
% input parameters
Sseo=0; Sno3o=0; Qo=0; Xso=0; Snh4o=0; Sno2o=0.2; So2o=0; KaL=0; Xaobo=0; Xnbo=0; Xnspo=0; Xhano=0; Xamxo=0;
%break
end
%withdrawal
%while t>t4 && t<=t5
if t>t4 && t<=t5
% input parameters
Sseo=0; Sno3o=0; Qo=600; Xso=0; Snh4o=0; Sno2o=0.2; So2o=0; KaL=0; Xaobo=0; Xnbo=0; Xnspo=0; Xhano=0; Xamxo=0;
%break
end
fval(1)=((Vo*Xaobo-Vo*Xaob/SRT)/V + u1*(So2/(Ko2aob+So2))*(Snh4/(Knh4aob+Snh4))*Xaob - baob*(So2/(Ko2aob+So2))*Xaob - baob*naob*((Ko2aob/(Ko2aob+So2))*(Sno2+Sno3)/(Kno3aob+Sno2+Sno3))*Xaob);
fval(2)=((Vo*Xnbo-Vo*Xnb/SRT)/V + u2*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xnb - bnb*(So2/(Ko2nb+So2))*Xnb - bnb*nnb*((Ko2nb/(Ko2nb+So2))*(Sno2+Sno3)/(Kno3nb+Sno2+Sno3))*Xnb);
fval(3)=((Vo*Xnspo-Vo*Xnsp/SRT)/V + u3*(Sno2/(Kno2nsp+Sno2))*(So2/(Ko2nsp+So2))*Xnsp - bnsp*(So2/(Ko2nsp+So2))*Xnsp - bnsp*nnsp*((Ko2nsp/(Ko2nsp+So2))*(Sno2+Sno3)/(Kno3nsp+Sno2+Sno3))*Xnsp);
fval(4)=((Vo*Xamxo-Vo*Xamx/SRT)/V + u5*((Snh4/(Knh4amx+Snh4))*(Sno2/(Kno2amx+Sno2))*(Ko2amx/(Ko2amx+So2)))*Xamx-bamx*(So2/(Ko2amx+So2))*Xamx - bamx*namx*((Ko2amx/(Ko2amx+So2))*(Sno2+Sno3)/(Kno3amx+Sno2+Sno3))*Xamx);
fval(5)=((Vo*Xhano-Vo*Xhan/SRT)/V + (u6*(Sse/(KSsehan+Sse))*(Sno2/(Kno2han+Sno2))*(Ko2han/(Ko2han+So2)) + u6*(Sse/(KSsehan+Sse))*(Sno3/(Kno3han+Sno3))*(Ko2han/(Ko2han+So2))+u6*((Sse/(KSsehan+Sse))*(So2/(Ko2han+So2)))*Xhan)- bhan*(So2/(Ko2han+So2))*Xhan - bhan*nhan*((Ko2han/(Ko2han+So2))*(Sno2+Sno3)/(Kno3han+Sno2+Sno3))*Xhan);
fval(6)=((Qo*Xso-Qo*Xs)/V - KH*((Xs/Xhan)/(KX+(Xs/Xhan)))*Xhan);
fval(7)=((Qo*Sseo-Qo*Sse)/V - (u6*nhan*(1/YNo2han)*(Sse/(KSsehan+Sse))*(Sno2/(Kno2han+Sno2))*(Ko2han/(Ko2han+So2)))*Xhan -(u6*nhan*(1/YNo3han)*(Sse/(KSsehan+Sse))*(Sno3/(Kno3han+Sno3))*(Ko2han/(Ko2han+So2)))*Xhan - ((u6/Yhaer)*(Sse/(KSsehan+Sse))*(So2/(Ko2han+So2))*Xhan) + KH*((Xs/Xhan)/(KX+(Xs/Xhan)))*Xhan);
fval(8)=((Qo*Sno3o-Qo*Sno3)/V - (u6*nhan*((1-YNo3han)/(2.86*YNo3han))*(Sse/(KSsehan+Sse))*(Sno3/(Kno3han+Sno3))*(Ko2han/(Ko2han+So2)))*Xhan + (u5/1.14)*(Snh4/(Knh4amx+Snh4))*(Sno2/(Kno2amx+Sno2))*(Ko2amx/(Ko2amx+So2))*Xamx + (u2/Ynb*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xnb) + ((u3/Ynsp)*(Sno2/(Kno2nsp+Sno2))*(So2/(Ko2nsp+So2))*Xnsp) - ((1-fI)/2.86)*baob*naob*(Ko2aob/(Ko2aob+So2))*(Sno2+Sno3)/(Kno3aob+Sno2+Sno3)*Xaob - ((1-fI)/2.86)*bnb*nnb*(Ko2nb/(Ko2nb+So2))*(Sno2+Sno3)/(Kno3nb+Sno2+Sno3)*Xnb - ((1-fI)/2.86)*bnsp*nnsp*(Ko2nsp/(Ko2nsp+So2))*(Sno2+Sno3)/(Kno3nsp+Sno2+Sno3)*Xnsp - ((1-fI)/2.86)*bamx*namx*(Ko2amx/(Ko2amx+So2))*(Sno2+Sno3)/(Kno3amx+Sno2+Sno3)*Xamx-((1-fI)/2.86)*bhan*nhan*(Ko2han/(Ko2han+So2))*(Sno2+Sno3)/(Kno3han+Sno2+Sno3)*Xhan);
fval(9)=((Qo*Sno2o-Qo*Sno2)/V - (((1-YNo2han)/(1.71*YNo2han))*u6*nhan*(Sse/(KSsehan+Sse))*(Sno2/(Kno2han+Sno2))*(Ko2han/(Ko2han+So2))*Xhan) - (((u5/Yamx)+(u5/1.14))*(Snh4/(Knh4amx+Snh4))*(Sno2/(Kno2amx+Sno2))*(Ko2amx/(Ko2amx+So2)))*Xamx + ((u1/Yaob)*(So2/(Ko2aob+So2))*(Snh4/(Knh4aob+Snh4)))*Xaob - ((u2/Ynb)*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xnb) - ((u3/Ynsp)*(Sno2/(Kno2nsp+Sno2))*(So2/(Ko2nsp+So2)))*Xnsp);
fval(10)=((Qo*Snh4o-Qo*Snh4)/V - (u1*(INBM+1/Yaob)*(So2/(Ko2aob+So2))*(Snh4/(Knh4aob+Snh4)))*Xaob -(u5*(INBM+1/Yamx)*(Snh4/(Knh4amx+Snh4))*(Sno2/(Kno2amx+Sno2))*(Ko2amx/(Ko2amx+So2)))*Xamx - (u2*INBM*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2)))*Xnb -(u3*INBM*(Sno2/(Kno2nsp+Sno2))*(So2/(Ko2nsp+So2)))*Xnsp - INBM*u6*nhan*((Sse/(KSsehan+Sse))*(Sno2/(Kno2han+Sno2))*(Ko2han/(Ko2han+So2)))*Xhan - u6*INBM*nhan*((Sse/(KSsehan+Sse))*(Sno3/(Kno3han+Sno3))*(Ko2han/(Ko2han+So2)))*Xhan - u6*INBM*((Sse/(KSsehan+Sse))*(So2/(Ko2han+So2)))*Xhan + (INBM-(fI*INXI))*baob*naob*((Ko2aob/(Ko2aob+So2))*(Sno2+Sno3)/(Kno3aob+Sno2+Sno3))*Xaob + (INBM-(fI*INXI))*bnb*nnb*((Ko2nb/(Ko2nb+So2))*(Sno2+Sno3)/(Kno3nsp+Sno2+Sno3))*Xnb +(INBM-(fI*INXI))*bnsp*nnsp*((Ko2nsp/(Ko2nsp+So2))*(Sno2+Sno3)/(Kno3nsp+Sno2+Sno3))*Xnsp + (INBM-(fI*INXI))*bamx*namx*((Ko2amx/(Ko2amx+So2))*(Sno2+Sno3)/(Kno3amx+Sno2+Sno3))*Xamx + (INBM-(fI*INXI))*bhan*nhan*((Ko2han/(Ko2han+So2))*(Sno2+Sno3)/(Kno3han+Sno2+Sno3))*Xhan +(INBM-(fI*INXI))*baob*(So2/(Ko2aob+So2))*Xaob + (INBM-(fI*INXI))*bnb*(So2/(Ko2nb+So2))*Xnb +(INBM-(fI*INXI))*bnsp*(So2/(Ko2nsp+So2))*Xnsp + (INBM-(fI*INXI))*bamx*(So2/(Ko2amx+So2))*Xamx + (INBM-(fI*INXI))*bhan*(So2/(Ko2han+So2))*Xhan);
fval(11)=((Qo*So2o-Qo*So2)/V + KaL*(8-So2) -(((1-Yhaer)/Yhaer)*u6*(Sse/(KSsehan+Sse))*(So2/(Ko2han+So2)))*Xhan -(((3.43-Yaob)/Yaob)*u1*(So2/(Ko2aob+So2))*(Snh4/(Knh4aob+Snh4)))*Xaob - (((1.14-Ynb)/Ynb)*u2*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2)))*Xnb - (((1.14-Ynsp)/Ynsp)*u3*(Sno2/(Kno2nsp+Sno2))*(So2/(Ko2nsp+So2)))*Xnsp - (1-fI)*baob*(So2/(Ko2aob+So2))*Xaob - (1-fI)*bnb*(So2/(Ko2nb+So2))*Xnb -(1-fI)*bnsp*(So2/(Ko2nsp+So2))*Xnsp - (1-fI)*bamx*(So2/(Ko2amx+So2))*Xamx - (1-fI)*bhan*(So2/(Ko2han+So2))*Xhan);
if fval(10)>maxSnh4
fval(10)=maxSnh4;
%else
% dy=fval(10);
end
dy = fval.';
end
%end
  2 Comments
KIPROTICH KOSGEY
KIPROTICH KOSGEY on 13 Jul 2023
Thank you for the insights. I will re-structure it further. However, I have coded for fval(10)<=200 but it is still going over this value as is evident in figure(1).
Please advise on this on this.
Torsten
Torsten on 13 Jul 2023
Edited: Torsten on 13 Jul 2023
The figure doesn't display "fval", but "y".
And the mass matrix must be of size 11x11, not 55x55.
But since you don't use the options settings, MATLAB does not throw an error.

Sign in to comment.

More Answers (0)

Categories

Find more on Labeling, Segmentation, and Detection in Help Center and File Exchange

Products


Release

R2023a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!