# How to control the outcome of an ODE?

4 views (last 30 days)
Danny Helwegen on 8 Mar 2021
Commented: Danny Helwegen on 8 Mar 2021
Hi everyone,
I have a system of ODE's that I solve by a ode15s solver. That all works fine, however, I want to try to controll the system. My current system is as follows:
Main script:
clear;clc;
%% Initial Conditions
p.CH_0 = 2000;
p.CO_0 = 1000;
p.T = 323.15;
p.Tres = 3140;
p.H = -400*1000;
%% Solver
tspan = [0 10000];
y0 = [p.CH_0;0;0;p.CO_0;p.T];
[times,ysolutions] = ode15s(@(t,y)func(t,y,p),tspan,y0);
plot(times,ysolutions(:,5));
xlabel('Time')
ylabel('Temperature [K]')
set(gca,'FontSize',12)
Ode's:
function dydt = func(t,y,p)
%% Allocate Memory
dydt = zeros(5,1);
%% Rates
Rate_1 = 16.6950 * 10e7 .* exp(-76*1000./(8.*y(5))).* y(1) .* y(4);
Rate_2 = 16.6950 * 70 .* exp(-46*1000./(8.*y(5))).* y(2) .* y(4);
%% Odes
dydt(1) = 2 .* ( (1./p.Tres) .* (p.CH_0 - y(1)) -Rate_1 );
dydt(2) = 2 .* ( -(1./p.Tres) .* y(2) + Rate_1 - Rate_2 );
dydt(3) = 2 .* ( -(1./p.Tres) .* y(3) + Rate_2 );
dydt(4) = 2.5 .* ( (10./p.Tres) .* (p.CO_0 - y(4)) - Rate_1 - Rate_2 );
% Temperature
dydt(5) = ( ( (1./p.Tres) * 2003360)*(p.T-y(5)) - p.H * (Rate_1 + Rate_2) ) / (1181545);
end
Giving the following figure for the temperature (dydt(5)): What I want to achieve by adding an additional term in the fifth ode is to controll the outlet temperature. I wan't to let the temperature rise up to 600 K, but not higher, I want to keep and controll it at this 600 K by varying the newly added term (U* (600 - y(5))). The rewritten fifth ode becomes:
dydt(5) = ( ( (1./p.Tres) * 2003360)*(p.T-y(5)) - p.H * (Rate_1 + Rate_2) + U * (600 - y(5)) ) / (1181545);
I have already tried to define U as an ODE, but this didn't work as I got some errors. Upon investigation I saw that it is perhaps possible to obtain the wanted result by making use of an event(?). But I don't know how to implement this.

Alan Stevens on 8 Mar 2021
How about just modifying dydt(5) within the function to
% Temperature
if y(5)>=600
dydt(5) = 0;
else
dydt(5) = ( ( (1./p.Tres) * 2003360)*(p.T-y(5)) - p.H * (Rate_1 + Rate_2) ) / (1181545);
end
Danny Helwegen on 8 Mar 2021
Ah yes, thats correct totally forgot about that part. As the specific heat removal mechanism isn't really important at the moment I think I can solve it by taking the difference between the steady state at 600 K and a reference temperature that isn't cooled down. Thanks for your help.