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.

조회 수: 2 (최근 30일)
KIPROTICH KOSGEY
KIPROTICH KOSGEY 2023년 7월 13일
편집: Torsten 2023년 7월 13일
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
이전 댓글 1개 표시
KIPROTICH KOSGEY
KIPROTICH KOSGEY 2023년 7월 13일
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 2023년 7월 13일
"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.

댓글을 달려면 로그인하십시오.

이 질문에 답변하려면 로그인하십시오.

채택된 답변

Torsten
Torsten 2023년 7월 13일
이동: Torsten 2023년 7월 13일
Ran in:
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
KIPROTICH KOSGEY
KIPROTICH KOSGEY 2023년 7월 13일
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 2023년 7월 13일
편집: Torsten 2023년 7월 13일
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.

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help Center 및 File Exchange에서 Programming에 대해 자세히 알아보기

태그

제품


릴리스

R2023a

Community Treasure Hunt

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

Start Hunting!

Translated by