How to control the outcome of an ODE?
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Danny Helwegen
le 8 Mar 2021
Commenté : Danny Helwegen
le 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 commentaires
Réponse acceptée
Alan Stevens
le 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 commentaires
Alan Stevens
le 8 Mar 2021
Modifié(e) : Alan Stevens
le 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.
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!