ODE triggered event does not stop

I have a modeled system with three limit situations, all three with two boundaries each (top and bottom). I have designed an integration function such that, should those boundaries be crossed, the integration stops and corrects itself. However, when obtaining value=0 inside the event function, the integration is not always stopped, and no value ie is obtained. Moreover, at times the value ie is both top and bottom boundaries (ie = [5;6] for example), which does not make sense.

function xx = integration(t_ode, data, cond0, options, fail_T, fail_valve)

t      = t_ode(1);
cont_t = 1;

ww_earth = [data.w_earth; data.w_earth; 0]; 

while t(end) < t_ode(end)
    
    [t, aux ,~ ,~ , ie] = ode15s(@ode_Gauss_J2_drag, t_ode(cont_t:end),...
                         cond0, options, data, ww_earth, fail_T, ...
                         fail_valve);
                      
    xx(cont_t:cont_t+length(t)-1,:) = aux;                                  % Saves ode solutions
    
    i1 = find(ie == 1 | ie == 2);
    i2 = find(ie == 3 | ie == 4);
    i3 = find(ie == 5 | ie == 6);
    
    if     i1 == 1
        xx(end,7) = data.g;
        xx(end,8) = 0;                                                      % Resets va condition to zero, accelerometer
        disp(ie)                                                            % has reached the maximum/minimum
    elseif i1 == 2                                                                 
        xx(end,7) = -data.g;
        xx(end,8) = 0;
        disp(ie)
    end
    
    if     i2 == 3
        xx(end,10) = 10*data.A0;
    elseif i2 == 4                                                          
        xx(end,10) = 0;                                                           % has reached the maximum
    end
    
    if     i3 == 5
        xx(end,11) = 0;                                                     % Resets v_valv condition to zero, valve has reached the maximum
    elseif i3 == 6
        xx(end,11) = 0;                                                     % Resets v_valv condition to zero, valve has reached the maximum
    end
    
    cond0  = xx(end,:);
    cont_t = cont_t + length(t);
    
end
end

The event function is as follows:

function [value, isterminal ,direction] = valve_event(~ ,xx ,data ,~ ,~ ,~)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
isterminal = [ 1 ; 1 ; 1 ;  1;  1; 1];
direction  = [-1 ; 1 ;-1 ;  1; -1; 1];

value      = [xx(7)  <  data.g
              xx(7)  > -data.g
              xx(10) <  data.A0*10
              xx(10) >  0
              xx(10) < data.A0*10 && xx(11) <= 0
              xx(10) > 0          && xx(11) >= 0];

end

The function integrated is the one below.

function [df, T] = ode_Gauss_J2_drag(~, xx, data, ww_earth, failure, failure_valv)
%% Variables

a     = xx(1);
e     = xx(2);
i     = xx(3);
% O     = xx(4);
o     = xx(5);
nu    = xx(6);
xa    = xx(7);
va    = xx(8);
Vout  = xx(9);
xvalv = xx(10);
vvalv = xx(11);
I     = xx(12);

if failure_valv == 1
    vvalv = 0;
end

%% Orbital and cartesian elements

[rr_sc, vv_sc] = kep2car(xx(1:6), data.muE);

rr = rr_sc';
vv = vv_sc';
r  = norm(rr);

%% Perturbing acceleration in cartesian coordinates (aap_car)

% J2 PERTURBATION

J2    = 0.00108263;   
Re    = data.Re;
ax_j2 = 3/2 * J2*data.muE*Re^2/r^4 * (rr(1)/r * (5*rr(3)^2/r^2 - 1));
ay_j2 = 3/2 * J2*data.muE*Re^2/r^4 * (rr(2)/r * (5*rr(3)^2/r^2 - 1));
az_j2 = 3/2 * J2*data.muE*Re^2/r^4 * (rr(3)/r * (5*rr(3)^2/r^2 - 3));
aa_j2 = [ax_j2; ay_j2; az_j2];                                              % Acceleration vector due to J2

% AIR DRAG PERTURBATION
r_earth = radius_earth(data.ra, data.rc, rr);
height  = r - r_earth;
rho     = atm_model(height);

vv_rel = vv;% - cross(ww_earth, rr);
v_rel  = norm(vv_rel);

A_sc       = data.A*1e-6;                                                   % Frontal area [km^2]
D          = 0.5*A_sc*data.cd*rho*v_rel^2;                                  % Drag force modulus [kg*km/s^2]
a_drag_sc  = -D/data.M;                                                     % Spacecraft's drag acceleration modulus [km/s]
aa_drag_sc = a_drag_sc * vv_rel/norm(vv_rel);                               % Acceleration vector due to air drag
aa_drag_sc = aa_drag_sc';

% BOUNDARY CONDITIONS
xacc_max  =  data.g;
xacc_min  = -data.g;
xvalv_max =  data.A0*10;
xvalv_min =  0;

% THRUST ACCELERATION

alpha  = 2*pi - 2*acos(1 - 2*xvalv/10/data.A0);                             % Angle defining geometry
A_valv = data.A0*(alpha - sin(alpha))/2/pi;                                 % Area varying as a circular conduct being opened
mdot   = data.cdis*A_valv*sqrt(2/data.R/data.T2*data.k/(data.k - 1))*...    % Mass flow rate (choked)
         data.p2*sqrt((2/(data.k + 1))^(2/(data.k - 1)) - (2/(data.k + ...
         1))^((data.k + 1)/(data.k - 1)));
     
if failure == 1
    T = 0;
else
    T = mdot*sqrt(2*data.e*data.DV/data.m_i);                               % Thrust generated by the ion thruster [N]
end

a_thrust  = T/data.M * 1e-3;                                                % [km/s^2]
aa_thrust = a_thrust * vv_rel/norm(vv_rel);
aa_thrust = aa_thrust';

% TOTAL PERTURBED ACCELERATION

aap_car_sc = aa_j2 + aa_drag_sc + aa_thrust;                                % In orbital cartesian coordinates

%% Components of the perturbing acceleration [at,an,ah]

t     = vv/norm(vv);
hh    = cross(rr, vv);
h0    = hh/norm(hh);
n0    = cross(h0, t);
A_mat = [t; n0; h0];

aap_tnh_sc = A_mat*aap_car_sc;                                              % Tangent, normal and angular momentum
at_sc      = aap_tnh_sc(1);
an_sc      = aap_tnh_sc(2);
ah_sc      = aap_tnh_sc(3);

%% Differential equations

% Useful parameters
b0     = a*sqrt(1 - e^2);
p0     = b0^2/a;
n0     = sqrt(data.muE/a^3);
h0     = n0*a*b0;
r0     = p0/(1 + e*cos(nu));
theta0 = nu + o;
v0     = sqrt(2*data.muE/r0 - data.muE/a);

% ORBITAL MECHANICS

df(1,1) = 2*a^2*v0/data.muE*at_sc;                                          % a
df(2,1) = 1/v0*(2*(e + cos(nu))*at_sc - r0/a*sin(nu)*an_sc);                % e
df(3,1) = r0*cos(theta0)/h0*ah_sc;                                          % i
df(4,1) = r0*sin(theta0)/(h0*sin(i))*ah_sc;                                 % O
df(5,1) = 1/(e*v0)*(2*sin(nu)*at_sc + (2*e + r0/a*cos(nu))*an_sc) - ...     % o
          r0*sin(theta0)*cos(i)/(h0*sin(i))*ah_sc;
df(6,1) = h0/r0^2 - 1/(e*v0)*(2*sin(nu)*at_sc + (2*e + r0/a*cos(nu))*...    % nu
          an_sc); 


% ACCELEROMETER

Vout_dot = -1/data.Cf*2*data.ep*data.Aa*(data.g^2 + xa^2)/(data.g^2 - ...   % d(Vout)/d(t)
            xa^2)^2*va*data.Vbias;

Vx = xa/data.g*data.Vbias;                                                  % Voltage generated by rotor
Vc = data.Kpa*Vout + data.Kda*Vout_dot;                                     % Control voltage     

DV1 = data.Vbias - Vc - 0.5*Vx;                                             % Increment of V at stator 1
DV2 = data.Vbias + Vc + 0.5*Vx;                                             % Increment of V at stator 2

Fe1 = 0.5*data.ep*data.Aa*DV1^2/(data.g - xa)^2;                            % Electric force exerted by stator 1
Fe2 = 0.5*data.ep*data.Aa*DV2^2/(data.g + xa)^2;                            % Electric force exerted by stator 2

D_acc = norm(a_drag_sc)*1e3*data.m;                                         % Drag exerted upon the accelerator [N]
T_acc = norm(a_thrust) *1e3*data.m;                                         % Thrust exerted upon the accelerator [N]

df(7,1) = va;                                                               % d(xa)/d(t)
df(8,1) = 1/data.m*(T_acc - D_acc - Fe1 + Fe2);                             % d(va)/d(t)
df(9,1) = Vout_dot;                                                         % d(Vout)/d(t)

% CONTROL VALVE

df(10,1) = vvalv;                                                           % d(xvalv)/d(t)
df(11,1) = 1/data.m_fcv*(data.Kfcv*(10*data.A0 - xvalv) - data.Ki*I - ...   % d(vvalv)/d(t)
            data.c*vvalv);
df(12,1) = data.Kpv*Vout_dot  + data.Kiv*Vout;                              % d(I)/d(t)

%failure mode: valve not working conditions
if failure_valv == 1
    df(10,1) = 0;
    df(11,1) = 0;
end

% if maximum valve stroke is reached, and velocity and acceleration are
% positive, acceleration must be zero.
if xvalv >= xvalv_max && vvalv>=0 && df(11)>0
    df(11) = 0;
end

if xvalv <= xvalv_min && vvalv<=0 && df(11)<0
    df(11) = 0;
end

end

Answer 1

[Steven Lord]

Your event function shouldn't use the relational operators itself. The only possible values your events can take on are 0 and 1 (false and true.) Instead your event function should return the distance between xx(7) and data.g (using the first event definition) and let the ODE solver determine when that distance (which can be positive, negative, or zero) crosses zero.

If you look at the ballode example, note that its event function doesn't return whether or not the ball is above the ground (true or false.) It returns the height above ground and the ODE solver determines when it hits and stops the solver at that point.

[Rafael Félix Soriano]

I have also tried this, but that would not let me use conditions 5 and 6 of the event. Instead, I've tried with if and changing value to 0 if the condition is met, but it does not work.

[Rafael Félix Soriano]

I don't know why, but for the same instant t, the event function does not use the same value of xx(10) than ode_Gauss_J2_drag does. Therefore, by the time valve_event detects the event, the actual value of xx(10) is already wrong (and has been for several iterations).

I've modified valve_event in several ways, but the error keeps appearing.

function [value, isterminal ,direction] = valve_event(t ,xx ,data ,~ ,~ ,~)

isterminal = [ 1;  1; 1;  1;  1;  1];
direction  = [ 1; -1; 1; -1; -1; -1];
% value      = [ 1;  1; 1;  1;  1;  1];

% if     xx(7) > data.g
%     value(1) = 0;
% elseif xx(7) < -data.g
%     value(2) = 0;
% end
% 
% if     xx(10) > 10*data.A0
%     value(3) = 0;
% elseif xx(10) < 0
%     value(4) = 0;
% end

value = [xx(7)  - data.g
         xx(7)  + data.g
         xx(10) - 10*data.A0
         xx(10)
         1
         1];

if     xx(10) >= 10*data.A0 && xx(11) > 0
    value(5) = 0;
elseif xx(10) <= 0          && xx(11) < 0
    value(6) = 0;
end

end

[Rafael Félix Soriano]

I am now trying with ode23s but, for some reason, the event is not checked at every iteration of the integrator. I have tried displaying the variable value inside valve_event and it does not appear at every iteration.