MATLAB Answers

How to control the outcome of an ODE?

4 views (last 30 days)
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.

Accepted Answer

Alan Stevens
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
Danny Helwegen
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.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by