How do I fix negative values in a function to zero ?
7 ビュー (過去 30 日間)
古いコメントを表示
I have the following script and function with a code to fix the negative values to zero but it is failing for some reason. What could be the issue?
%Script:
%initial conditions
y0=[18.16;0.52;25.64;172.6;179.48;0.0005;2;0.0005;55;50;0.5];
tSpan=[0 535];
M1=eye(11);
options = odeset('RelTol',1e-20,'AbsTol',1e-10,'Mass',M1,'MassSingular', 'maybe','MStateDependence','strong','InitialStep',1e-10, 'MStateDependence','strong');
[tSol, ySol]=ode15s(@(t,y) MBBR12(t,y), tSpan, y0, options);
%Funtion as follows:
function fval=MBBR12(t,y)
%Define the three variables
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.047*exp(O3*(T+293)); Ko2nsp=2.2; Kno2nsp=5.5; Ynsp=0.041;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=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=75; Xso=250;
Snh4o=650; Sno2o=0.2; So2o=0;
%oxygen transfer parameters
KaL=80;
V=1400; SRT=40;
%fix output values to zero
Sno3(Sno3<0)=0;Xaob(Xaob<0)=0;Xnb(Xnb<0)=0;Xnsp(Xnsp<0)=0;Xamx(Xamx<0)=0;Xhan(Xhan<0)=0;Xs(Xs<0)=0;
Sse(Sse<0)=0;Sno2(Sno2<0)=0;Snh4(Snh4<0)=0;So2(So2<0)=0;
% main functions CSTR 1
fval=[((0-V*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)
((0-V*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)
((0-V*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)
((0-V*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)
((0-V*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)
((Qo*Xso-Qo*Xs)/V - KH*((Xs/Xhan)/(KX+(Xs/Xhan)))*Xhan)
((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)
((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)
((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)
((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)
((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)];
end
0 件のコメント
採用された回答
John D'Errico
2023 年 7 月 9 日
Have something like this be the last thing you do in the function.
x = max(x,0);
Or,
x(x<0) = 0;
3 件のコメント
Steven Lord
2023 年 9 月 1 日
For this particular scenario, where you're trying to constrain the solution to a system of ODEs to be non-negative, rather than calling max at the end of the ODE function I'd use the NonNegative option for ODEs if it's supported for the ODE solver you're using. See this example in the documentation.
その他の回答 (0 件)
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!