Event function with DAE system

2023 10 2
1 个回答
回答已采纳
更新时间：2023 10 27
5 次查看（30 天）

0 个投票

I have a DAE system for solving a set of diffrential and algebraic equations(4 differential equations and 13 algebraic). I am using the ode15s solver. In my system, I am calculating the Pressure of a vapour in a storage tank which is increasing with time at the moment.
Once the pressure reaches an upper threshold value, I want 2 of the 4 differential equations to change slightly. These updated equations will make the pressure drop and once the pressure reaches the lower threshold value, the system should revert back to the orginal differential equations.
Any idea how I can do this using event functions??
This is how I am defining my DAE_system function with a mass matrix
options = odeset('Mass', M, 'MassSingular', 'yes', 'RelTol', 1e-4, 'AbsTol', 1e-6);
[t, y] = ode15s(@(t,y) DAE_System(t,y,Params,userInput_geo), tspan, y0, options);
These are the original equations:
dy_dt(2,1) = Q_av - Q_v + M_e*hv + Pv*dy_dt(1, 1);
dy_dt(3,1) = M_e;
These are the equations I want once the upper threshold limit of Pv is reached:
dy_dt(2,1) = Q_av - Q_v + M_e*hv + Pv*dy_dt(1,1) - M_va*hv;
dy_dt(3,1) = M_e - M_va;

请先登录，再回答此问题。

 采纳的回答

Torsten
Torsten 2023-10-2
编辑：Torsten 2023-10-2

1 个投票

options = odeset('Mass', M, 'MassSingular', 'yes', 'RelTol', 1e-4, 'AbsTol', 1e-6,'Events',@event);
function [position,isterminal,direction] = event(t,y)
Pvmax = ...;
position = y(?) - Pvmax; % The value that we want to be zero, set the ? to the index of y that contains Pv
isterminal = 1; % Halt integration
direction = 0; % The zero can be approached from either direction
end

17 个评论
显示 15更早的评论 隐藏 15更早的评论

Yash Vardhan
Yash Vardhan 2023-10-2
Makes sense: But how do I ensure the solver uses the new equations and allow it to revert back to the original equations once the pressure drops?
Torsten
Torsten 2023-10-2
Integrate with the original equations until the event happens. Use the results from the first run as new initial conditions for a call of ode15s with the new equations.
I hope your system doesn't oscillate between Pv < Pvmax and Pv > Pvmax several times.
Yash Vardhan
Yash Vardhan 2023-10-2
@Torsten Are you suggesting I make a new DAE_system function with the new equations and use the last values of the state variables from the first event occurance as the initial conditions for the ew DAE_system ??
Torsten
Torsten 2023-10-2
编辑：Torsten 2023-10-2
Exactly. You have to use two DAE_system functions. Or you use one function and pass a flag to the function in order to indicate which phase of the calculation you are in.
Something like
iflag = 1;
options1 = odeset('Mass', M, 'MassSingular', 'yes', 'RelTol', 1e-4, 'AbsTol', 1e-6,'Events',@event);
[t1, y1] = ode15s(@(t,y) DAE_System(t,y,Params,userInput_geo,iflag), tspan, y0, options1);
iflag = 2;
y0 = y1(end,:);
tspan = [t1(end) tspan(end)];
options2 = odeset('Mass', M, 'MassSingular', 'yes', 'RelTol', 1e-4, 'AbsTol', 1e-6);
[t2, y2] = ode15s(@(t,y) DAE_System(t,y,Params,userInput_geo,iflag), tspan, y0, options2);
t = [t1;t2(2:end)];
y = [y1;y2(2:end,:)];
plot(t,y)
function dydt = DAE_System(t,y,Params,userInput_geo,iflag)
...
if iflag == 1
dy_dt(2,1) = Q_av - Q_v + M_e*hv + Pv*dy_dt(1, 1);
dy_dt(3,1) = M_e;
elseif iflag == 2
dy_dt(2,1) = Q_av - Q_v + M_e*hv + Pv*dy_dt(1,1) - M_va*hv;
dy_dt(3,1) = M_e - M_va;
end
...
end
Yash Vardhan
Yash Vardhan 2023-10-3
@Torsten Thanks! That seems to work but any idea how can I make the system revert back to the original equations if Pv drops below Pv_min?
Torsten
Torsten 2023-10-3
Use a second event function (now with Pv_min instead of Pv_max) and restart ...
Yash Vardhan
Yash Vardhan 2023-10-3
@Torsten Actually I dont want the solver to restart but to just start using the original equations again from the time when Pv <= Pvmin
Torsten
Torsten 2023-10-3
编辑：Torsten 2023-10-3
Yes, that's what I meant with "restart".
t = [];
y = [];
tend = -1;
while tend < tspan(end)
iflag = 1;
options1 = odeset('Mass', M, 'MassSingular', 'yes', 'RelTol', 1e-4, 'AbsTol', 1e-6,'Events',@event_Pvmax);
[t1, y1] = ode15s(@(t,y) DAE_System(t,y,Params,userInput_geo,iflag), tspan, y0, options1);
iflag = 2;
y0 = y1(end,:);
tspan = [t1(end) tspan(end)];
options2 = odeset('Mass', M, 'MassSingular', 'yes', 'RelTol', 1e-4, 'AbsTol', 1e-6,'Events',@event_Pvmin);
[t2, y2] = ode15s(@(t,y) DAE_System(t,y,Params,userInput_geo,iflag), tspan, y0, options2);
y0 = y2(end,:);
tspan = [t2(end) tspan(end)];
t = [t;t1;t2(2:end)];
y = [y;y1;y2(2:end,:)];
tend = t(end);
end
Yash Vardhan
Yash Vardhan 2023-10-5
@Torsten
Thanks a lot! That worked wonders!
Yash Vardhan
Yash Vardhan 2023-10-5
@Torsten
If there are more parameters I want to have event functions for. Is there a way of combining the event functions?
For eg: I have a variable M_l and I dont want it to drop to 0 and also not exceed M_lmax as well.
Torsten
Torsten 2023-10-5
编辑：Torsten 2023-10-5
You can, and the parameter "ie" in the list of output parameters from the ODE solver tells you which event has occured:
[t,y,te,ye,ie] = odeXY(odefun,tspan,y0,options)
For more information, consult
https://uk.mathworks.com/help/matlab/math/ode-event-location.html
Torsten
Torsten 2023-10-9
I'd try
direction = [0 0]
And the output list must read
[t_event, y_event,~,~,ie] = ode15s(@(t, y) DAE_System(t, y, Params, userInput_geo, event_flag), tspan, y0, options);
instead of
[t_event, y_event,ie, ~, ~] = ode15s(@(t, y) DAE_System(t, y, Params, userInput_geo, event_flag), tspan, y0, options);
If this doesn't work, it's necessary to have executable code.
Yash Vardhan
Yash Vardhan 2023-10-27
编辑：Yash Vardhan 2023-10-27
@Torsten
I am now using the ode15i solver and am running into an issue with the tspan and the y [ ] matrix.
I have 18 state variables and the initial values are defined in y0 and the yp0 has the derivatives.
y0 = [......];
yp0 = y0*0;
First of all the ie value I am getting is an array with random values (eg: 50 10 241 3 33 102) and not a single number. With the ode15s solver, I was getting a single value based on which error was occuring.
Secondly, the solver runs in a very weird manner. It goes to the the Pv_max at a certain time t1. Then instead of charging the equations it uses, it starts the simulation from scratch with the inital values and t0. This keeps on happening again and again.
This is my while loop for running the solver:
%% Running the Solver with event functions
event_flag = 1;
[y0,yp0] = decic(@(t,y,yp)DAE_System(t,y,yp,Params,userInput_geo,event_flag),0,y0,[1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0],yp0,[]);
while true
options = odeset('RelTol', 1e-2, 'Events', @(t, y, yp) event_function2(t, y, T_crit, P_crit, V_T, Pv_ul, Pv_ll));
[t_event, y_event, yp_event, te, ye, ie] = ode15i(@(t, y,yp) DAE_System(t, y, yp, Params, userInput_geo, event_flag), tspan, y0,yp0,options);
% Making 1 final array for t and 1 array for y and concatenating the results
t = [t; t_event];
y = [y; y_event];
if isempty(ie) % No event occurred, break the loop
break;
elseif any(ie == [1,2,3,4,5]) % End simulation
break;
end
elseif any(ie == [6,7])
if event_flag == 1
event_flag = 2;
else
event_flag = 1;
end
% Updating tspan, y0 and yp0 for continuing the simulation
if t_event(end) < tspan_finish
tspan = [t_event(end), tspan_finish];
y0 = y_event(end, :);
yp0 = yp_event(end, :);
else
break;
end
% Printing which event occured when
fprintf('Event occurred at t = %f. event_flag set to %d\n', t_event(end), event_flag);
end
end
This is my event function:
%% Event Function
function [check,isterminal,direction] = event_function2(t, y, yp, T_crit, P_crit, V_T, Pv_ul, Pv_ll)
condition1 = y(6);
condition2 = y(1) - V_T*0.999;
condition3 = y(7) - T_crit;
condition4 = y(8) - T_crit;
condition5 = y(9) - P_crit;
condition6 = y(9) - Pv_ul;
condition7 = Pv_ll - y(9);
check = [condition1; condition2; condition3; condition4; condition5; condition6; condition7];
isterminal = [1;1;1;1;1;1;1]; % terminating intergartion when event occurs
direction = [0;0;0;0;0;1;1]; % use default settings for direction
end
Could you help me get this event function and event handling to work??
Thanks!
Torsten
Torsten 2023-10-27
if isempty(ie) % No event occurred, break the loop
break;
elseif any(ie == [1,2,3,4,5]) % End simulation
break;
end
elseif any(ie == [6,7])
Doesn't this if-construct throw an error ?
Yash Vardhan
Yash Vardhan 2023-10-27
@Torsten
Yes It was throwing an error which is why I changed the system a little bit
while true
% Changing the options to use the correct event_function based on event_flag
[y0,yp0] = decic(@(t,y,yp)DAE_System(t,y,yp,Params,userInput_geo,event_flag),0,y0,[1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0],yp0,[]);
options = odeset('Events', @(t, y, yp) event_function2(t, y, yp, T_crit, P_crit, V_T, Pv_ul, Pv_ll));
[t_event, y_event, ~, ~, ie] = ode15i(@(t, y,yp) DAE_System(t, y, yp, Params, userInput_geo, event_flag), tspan, y0,yp0,options);
% Making 1 final array for t and 1 array for y and concatenating the
% results after every event
t = [t; t_event];
y = [y; y_event];
if isempty(ie) % No event occurred, break the loop
break;
elseif ie == 1||ie == 2 || ie == 3||ie == 4 || ie == 5 % End simulation
break;
elseif ie == 6
if event_flag == 1
ie
event_flag = 2;
end
% Updating tspan and y0 for continuing the simulation
if t_event(end) < tspan_finish
tspan = [t_event(end), tspan_finish];
y0 = y_event(end, :);
else
break;
end
% Printing which event occured when
fprintf('Event occurred at t = %f. event_flag set to %d\n', t_event(end), event_flag);
elseif ie == 7
ie
event_flag = 1;
if t_event(end) < tspan_finish
tspan = [t_event(end), tspan_finish];
y0 = y_event(end, :);
else
break;
end
% Printing which event occured when
fprintf('Event occurred at t = %f. event_flag set to %d\n', t_event(end), event_flag);
end
end
Now the the issue I am getting is with regards the tolerances
Warning: Failure at t=2.649540e+04. Unable to meet integration tolerances without reducing the step size below the
smallest value allowed (9.413057e-11) at time t.
> In ode15i (line 398)
Any clue?
Torsten
Torsten 2023-10-27
You don't supply the actual derivatives to "decic" when you make a restart.
The yp0 in the call
decic(@(t,y,yp)DAE_System(t,y,yp,Params,userInput_geo,event_flag),0,y0,[1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0],yp0,[]);
is still the yp0 from the last t_event(end), not from the actual one.
After an event happened, use "deval" to recompute yp0 from t_event(end) and y_event(end,:).
[~,yp] = deval( )
https://uk.mathworks.com/help/matlab/ref/deval.html
Further, I'm not sure whether
elseif ie == 7
ie
event_flag = 1;
suffices or something as in the case ie = 6:
if event_flag == ...
ie
event_flag = ...;
end
The error message indicates that ode15i couldn't continue integration because it could no longer retain the error tolerances. Reasons can be a singularity in the solution, a programming error, too strict or weak tolerances,... . The only way to cope with it is to look at the solution at t=2.649540e+04 or slightly before this time and try to find the reason for the failure of the integrator.

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息

产品

版本

R2022b

标签

Community Treasure Hunt

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

Start Hunting!

Translated by