Event function with DAE system

Question

Yash Vardhan 2023-10-2

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2028219-event-function-with-dae-system

评论： Torsten 2023-10-27

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;

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Torsten 2023-10-2

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2028219-event-function-with-dae-system#answer_1323584

编辑：Torsten 2023-10-2

在 MATLAB Online 中打开

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 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

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!

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

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 2023-10-27

在 MATLAB Online 中打开

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.

请先登录，再进行评论。

Event function with DAE system

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

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

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Event function with DAE system

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

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

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

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