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).

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

Answer 1

[Walter Roberson]

You already have

[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.

[Tom Goodland]

Thanks all the graphs load for the different temperatures, but they all load very quickly lasting a few seconds as every split second a new graph shows, but when all have been loaded the rest dissapear and i'm left with 3 graphs for the original 455K problem. Do you know a way to make it so all the graphs stay instead of deleting themselves once loaded?

[Tom Goodland]

Is there something I can do with the plot function and the thermostat_setpoints?

[Walter Roberson]

You can use hold on

msg = compose("@%.2f", setpoint);
figure(1)
H = plot(t, [ca,cb,cs,cd]);
H(1).DisplayName = "ca" + msg;
H(2).DisplayName = "cb" + msg;
H(3).DisplayName = "cc" + msg;
H(4).DisplayName = "cd" + msg;
hold on
legend show

figure(2)
plot(t, P, 'displayname', "P" + msg);
hold on
legend show

figure(3)
plot(t, T_2, 'displayname', "T_2" + msg);
legend show

end

[Tom Goodland]

When I run it matlab says Unrecognized function or variable 'setpoint'. Error in twoc (line 19)

msg = compose("@%.2f", setpoint). Do you know how to fix this?

[Image Analyst]

Maybe it was supposed to be current_setpoint instead.

[Stephen23]

"what is the @%,2f part?"

Using a dot (not a comma) specifies the number of decimal digits after the radix symbol: https://www.mathworks.com/help/matlab/ref/compose.html#bu6eywi_sep_bu6eywi-formatSpec

compose("@%.2f", [2/3,pi,1234567890])

ans = 1×3 string array
    "@0.67"    "@3.14"    "@1234567890.00"

[Tom Goodland]

Do you know how I can write some code to stop the for loop from doing its job of increasing the temperature then running the ode again when T_2 = 500K? I have tried with a break but I keep on getting errors.

function twoc
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
thermostat_setpoints = 455:2.5:500;
ntemps = length(thermostat_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);
    t = [t1;te;t2];
    y = [y1;ye;y2];

[Tom Goodland]

@Walter Roberson@Thorsten

Do you know how I can write some code to stop the for loop from doing its job of increasing the temperature then running the ode again when T_2 = 500K for any of the iterations? I have tried with a break but I keep on getting errors. I'd really appreciate any help.

[Walter Roberson]

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