How to stop ODE fuction in a certain value

6 次查看（过去 30 天）
Berk Han
Berk Han 2022-6-17
编辑： Torsten 2022-6-17
This is a basically batch bioreactor but I want a fed-batch reactor. In the code, I want to stop the ODE solver when my y(2) fuction values equal 20 or less than 20. Then, I want to continoue to solve ODE with again setting y(2) as its inial point y2o(=31). So basically I want to go back to inital data when ever y(2)=<20
Is there anyway to add break if else or for loop to ode45
function [tsoln, ysoln] = odemultiple
kd=0.076
Bmax=15.658
ib=0.616
K1=171.492
Umaxg=0.730
Ksg=1.293
Ksx=4.469
Umaxx=0.615
Yxsg=0.523
Yxsx=0.058
Ybsg=1.196
Ybsx=0.232
Ybxg=Ybsg/Yxsg
Ybxx=Ybsx/Yxsx
to = 0; tn = 100;
y1o = 0.9; y2o = 31; y3o = 32; y4o=0;
[tsoln,ysoln] = ode45(@odefun,[to tn],[y1o y2o y3o y4o]);
plot(tsoln,ysoln)
legend('y1','y2','y3','y4');
function dy = odefun(t,y)
Usg=(Umaxg*y(2))/((Ksg+y(2))*((1+y(3))/Ksx))
Usx=(Umaxx*y(3))/((Ksx+y(3))*((1+y(2))/Ksg))
Ug=(Usg+Usx)*(K1/(K1+y(2)+y(3)))*(1-(y(4)/Bmax))^(ib)
Unet=Ug-kd
dy = zeros(4,1);
dy(1) = Unet*y(1);
dy(2) = -Usg*((1/Yxsg)+(1/Ybsg))*y(1);
dy(3) = -Usx*((1/Yxsx)+(1/Ybsx))*y(1);
dy(4)=(Usg*Ybxg+Usx*Ybxx)*y(1);
end
end

请先登录，再回答此问题。

采纳的回答

patrick1704
patrick1704 2022-6-17
编辑：patrick1704 2022-6-17
From my perspective the already given answer only solves the problem of stopping the simulation. However, it does not restart it as also desired (if I understand the question correctly).
From my perspective, you may only solve this via "recursion" in a multiple-shooting kind of fashion. In an (inefficient) implementation, this would look like the following:
storeTime = [];
storeSol = [];
while true
[tsoln,ysoln] = ode45(@odefun,[to tn],[y1o y2o y3o y4o]);
if all(ysoln(:,2) > 20)
storeTime = [storeTime; tsoln];
storeSol = [storeSol; ysoln];
break;
else
idx = find(real(ysoln(:,2)) <= 20, 1, 'first');
to = tsoln(idx);
y1o = real(ysoln(idx, 1));
y3o = real(ysoln(idx, 3));
y4o = real(ysoln(idx, 4));
storeTime = [storeTime; tsoln(1:1:idx-1)];
storeSol = [storeSol; ysoln(1:1:idx-1,:)];
end
end
plot(storeTime,storeSol)
legend('y1','y2','y3','y4');
Naturally, you can also add the "Events" function to stop the simulation and reduce your computational effort.
This creates the desired restart of the simulation until the whole time history is above the desired threshold.
  2 个评论
Berk Han
Berk Han 2022-6-17
This is excatly what I want to do. Thank you so much.
I have one more question. If I want to loob this should I just add an for loop on top of this code?
patrick1704
patrick1704 2022-6-17
No problem.
I don't know what you want to do exactly but you can naturally put a for-loop "around" the while loop to e.g., evaluate it for different initial conditions.

请先登录，再进行评论。

更多回答（1 个）

Torsten
Torsten 2022-6-17
编辑：Torsten 2022-6-17
Ran in:
Better you use the Event facility of the ODE integrators.
Note that your results become complex-valued because y(4) becomes greater than Bmax. So y(4) has also be adapted for the integration beginning where y(2) <= 20.
odesolve()
Warning: Imaginary parts of complex X and/or Y arguments ignored.
function odesolve
kd=0.076;
Bmax=15.658;
ib=0.616;
K1=171.492;
Umaxg=0.730;
Ksg=1.293;
Ksx=4.469;
Umaxx=0.615;
Yxsg=0.523;
Yxsx=0.058;
Ybsg=1.196;
Ybsx=0.232;
Ybxg=Ybsg/Yxsg;
Ybxx=Ybsx/Yxsx;
to = 0; tn = 100;
y1o = 0.9; y2o = 31; y3o = 32; y4o=0;
options = odeset('Events',@restart);
[tsoln1,ysoln1,te,ye,ie] = ode45(@odefun,[to tn],[y1o y2o y3o y4o],options);
[tsoln2,ysoln2] = ode45(@odefun,[te tn],[ye(1) y2o ye(3) ye(4)]);
plot([tsoln1;tsoln2],[ysoln1;ysoln2])
legend('y1','y2','y3','y4')
function dy = odefun(t,y)
Usg=(Umaxg*y(2))/((Ksg+y(2))*((1+y(3))/Ksx));
Usx=(Umaxx*y(3))/((Ksx+y(3))*((1+y(2))/Ksg));
Ug=(Usg+Usx)*(K1/(K1+y(2)+y(3)))*(1-(y(4)/Bmax))^(ib);
Unet=Ug-kd;
dy = zeros(4,1);
dy(1) = Unet*y(1);
dy(2) = -Usg*((1/Yxsg)+(1/Ybsg))*y(1);
dy(3) = -Usx*((1/Yxsx)+(1/Ybsx))*y(1);
dy(4)=(Usg*Ybxg+Usx*Ybxx)*y(1);
end
function [value,isterminal,direction] = restart(t,y)
value = y(2)-20.0;
isterminal = 1;
direction = 0;
end
end

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by