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.
1 view (last 30 days)
Show older comments
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
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.
Accepted Answer
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);
%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
More Answers (0)
See Also
Categories
Find more on Labeling, Segmentation, and Detection in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!