How to control the outcome of an ODE?

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.

Thanks in advance.

Answer 1

[Alan Stevens]

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]

First of all, thanks for your answer which works. However, I also want to know what the value of U is overtime. So how much heat should I take away to stay at this temperature over time.

[Alan Stevens]

What's your physical heat removal mechanism? If you just put U*(600 - y(5)) then this will be zero if y(5) stays at 600. You need something like U*(y(5) - Tsink), perhaps.

[Danny Helwegen]

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.