ODE triggered event does not stop

5 次查看（过去 30 天）
Rafael Félix Soriano
编辑： Rafael Félix Soriano 2019-12-13
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

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

回答（1 个）

Steven Lord
Steven Lord 2019-12-11
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.
  3 个评论
显示 1更早的评论
Rafael Félix Soriano
Rafael Félix Soriano 2019-12-12
编辑：Rafael Félix Soriano 2019-12-13
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
Rafael Félix Soriano 2019-12-13
编辑：Rafael Félix Soriano 2019-12-13
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息

标签

产品


版本

R2019b

Community Treasure Hunt

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

Start Hunting!

Translated by