Event function with DAE system

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;

Answer 1

[Torsten]

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

[Yash Vardhan]

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]

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]

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

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]

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

Use a second event function (now with Pv_min instead of Pv_max) and restart ...

[Yash Vardhan]

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

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]

@Torsten

Thanks a lot! That worked wonders!

[Yash Vardhan]

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

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]

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]

@Torsten

Okay Let me try this!

[Yash Vardhan]

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

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]

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

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.