How to control the outcome of an ODE?
1 view (last 30 days)
Show older comments
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.
Thanks in advance.
0 Comments
Accepted Answer
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
3 Comments
Alan Stevens
on 8 Mar 2021
Edited: Alan Stevens
on 8 Mar 2021
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.
More Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!