# How to run the same code over and over with one different detail of an event changing each time (temperature the cooling water turns on).

3 views (last 30 days)
Tom Goodland on 17 Jan 2022
Commented: Walter Roberson on 19 Jan 2022
I have to find out the maximum time the cooling water can be turned off before being reinstated to avoid the reactor exploding (happens at T = 500 K or P = 45 atm). I want to use the backbone of the code I used for the previous part of my work and have a function keep on picking temperatures to turn on the cooling water above 455 K. The objective for me is to find out how late I can turn the cooling water on and the reactor not explode. Does anyone know what function or code I could use to have UA turn from 0 (NO COOLING WATER) to 2.77E6 (COOLING WATER ON) running the code again for each different temperature the water is turned on? (represents cooling water in the dT derivative)
I think a loop of some sort would help me run this code over and over again with the UA part of the dT/dt derivative changing at different temperatures. However, i'm very unfamiliar with loops and I would appreciate any help.If you need me to provide the code I've been using for the ode, the derivatives, variables, constants etc. let me know. I've posted below the code I used to make an event for previous work on the same problem where the cooling water turns on at 455K.
Any help would be greatly appreciated, thanks.
function main
tspan = [0 10]; %just as an example
initial_conditions = [4.3 5.1 3 0 4.4 422]; %its the values of the variables you are solving for at t=0
options = odeset('Events',@myEventsFcn);
iflag = 1;
[t1,y1,te,ye,ie] = ode15s(@(t,y)fun(t,y(1),y(2),y(3),y(4),y(5),y(6),iflag), tspan, initial_conditions,options);
iflag = 2;
[t2,y2] = ode15s(@(t,y)fun(t,y(1),y(2),y(3),y(4),y(5),y(6),iflag), [te tspan(end)], ye);
t = [t1;te;t2];
y = [y1;ye;y2];
%output t is your time range, and y an array with all your variables where:
ca = y(:,1);
cb = y(:,2);
cs = y(:,3);
cd = y(:,4);
P = y(:,5);
T_2 = y(:,6);
figure(1)
plot(t,[ca,cb,cs,cd])
figure(2)
plot(t,P)
figure(3)
plot(t,T_2)
end
function [value,isterminal,direction] = myEventsFcn(t,y)
value = y(6) - 455; % The value that we want to be zero
isterminal = 1; % Halt integration
direction = 0; % The zero can be approached from either direction

Walter Roberson on 17 Jan 2022
[t1,y1,te,ye,ie] = ode15s(@(t,y)fun(t,y(1),y(2),y(3),y(4),y(5),y(6),iflag), tspan, initial_conditions,options);
so you already know how to construct anonymous functions to pass in additional parameters.
You can do the same thing for the target temperature
thermostat_setpoints = 455:2.5:500;
ntemps = length(termostat_setpoints);
for tempidx = 1 : ntemps
current_setpoint = thermostat_setpoints(tempidx);
options = odeset('Events', @(t,y)myEventsFcn(t,y,current_setpoint));
[t1,y1,te,ye,ie] = ode15s(@(t,y)fun(t,y(1),y(2),y(3),y(4),y(5),y(6),1), tspan, initial_conditions, options);
[t2,y2] = ode15s(@(t,y)fun(t,y(1),y(2),y(3),y(4),y(5),y(6),2), [te tspan(end)], ye);
%and so on
end
function [value,isterminal,direction] = myEventsFcn(t,y,setpoint)
value = y(6) - setpoint; % The value that we want to be zero
isterminal = 1; % Halt integration
direction = 0; % The zero can be approached from either direction
end
I would further point out that if you were to add an additional event that detected "the reactor exploding", especially in combination with the second ode15s call, then instead of doing a for loop, you could do a binary search.
Walter Roberson on 19 Jan 2022
stop_temperature = 500;
thermostat_setpoints = 455:2.5:stop_temperature*(1-eps);
ntemps = length(termostat_setpoints);
for tempidx = 1 : ntemps
current_setpoint = thermostat_setpoints(tempidx);
options = odeset('Events', @(t,y)myEventsFcn(t,y,[stop_temperature,current_setpoint);
[t1,y1,te,ye,ie] = ode15s(@(t,y)fun(t,y(1),y(2),y(3),y(4),y(5),y(6),1), tspan, initial_conditions, options);
if any(ie == 1)
%index 1 --> explosion
break;
end
[t2,y2,te2,ye2,ie2] = ode15s(@(t,y)fun(t,y(1),y(2),y(3),y(4),y(5),y(6),2), [te tspan(end)], y1(end,:), options);
if any(ie2 == 1)
%index 1 --> explosion
break;
end
end
function [value,isterminal,direction] = myEventsFcn(t,y,setpoints)
value = y(6) - setpoints; % The value that we want to be zero
isterminal = [1 1]; % Halt integration
direction = [0 +1]; %expode --> always halt; turn on cooling only if temperature is increasing
end

R2021b

### Community Treasure Hunt

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

Start Hunting!

Translated by