How to stop ODE fuction in a certain value

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

Answer 1

[patrick1704]

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.

[Berk Han]

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]

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.

Answer 2

[Torsten]

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