I am having trouble triggering an event

2 次查看（过去 30 天）
Dean
Dean 2013-6-17
I want to trigger an event to signal a takeoff. The code plots the inrun for a snowboard jump and I am using an event to stop the ode when the "rider" takesoff. The code is posted below as well as the inrun and event files that the ode15 calls. The problem seems like the trigger extends past the data into NaN territory and then doesn't know what to do. For now, I have extended the data by 10 meters and triggered the event at the last point minus 10 meters (which is the real takeoff point). How do I make this work without temporarily extending the data?
% File calculates the rider position along the x axis and the velocity along
% the slope(i.e. sdot) and the acceleration along the slope(sdoubledot) for
% a time interval tspan.
%
% Ap = approach
% Tr = transition
% To = takeoff
% LengthAp = approach slope length
% LengthTo = takoff slope length
% R = radius of transition
% Inr = inrun
% Ps = parent slope
% imp = impact
close all;
clear all;
clc;
%format long
format short
global xInrun yInrun xPs yPs dydxInrun kurvature xExtInrun yExtInrun ext
global m CdA mu eta g h
global xland yland xLanding yLanding
global thetaAp thetaTo thetaPs
global xImpact yImpact vImpact
global xInr yInr
ft2m=12/39.37; % feet to meter conversion factor
m2ft=1/ft2m; % meters to feet conversion factor
deg2rad=pi/180; % degree to radian conversion factor
mps2mph=2.237; % meters per second to miles per hour conversion
mph2mps=1/mps2mph; % miles per hour to meters per second conversion
% given variables/entered variables
CdA = 0.279 % m^2 0.279 is least drag position
rho0 = 1.1839;
T0=298.15;
Y0=8772.;
T=0+T0;
Y=3870; % elevation of Loveland in meters
rho=rho0*exp(-T0*Y/T/Y0)*T0/T
g = 9.81
m = 60
mu = 0.1
eta = (CdA*rho)/(2*m)
thetaPs = 21*deg2rad;
thetaAp = thetaPs;
thetaTo = 20.6*deg2rad;
xRiderStart = 50
h=1.125;
xPs = 0:1:500;
yPs = -xPs*tan(thetaPs);
%%Calculating (x,y) from rider starting position to takeoff
NUMPOINTS = 100;
LengthAp = 60; % length of the approach-from start to transition
LengthTo = 3; % length of the takeoff
RadiusTr = 9; % radius of the transition segment
yRiderStart = interp1(xPs, yPs, xRiderStart); % interpolate y position of rider start
xTrStart = xRiderStart + LengthAp*cos(thetaAp);
yTrStart = yRiderStart-LengthAp*sin(thetaAp);
xToStart = xTrStart + RadiusTr*sin(thetaAp) + RadiusTr*sin(thetaTo);
yToStart = yTrStart - (RadiusTr*cos(thetaTo) - RadiusTr*cos(thetaAp));
xToEnd = xToStart + LengthTo*cos(thetaTo);
yToEnd = yToStart + LengthTo*sin(thetaTo);
% (x,y) data for approach
xAp = linspace(xRiderStart,xTrStart,NUMPOINTS);
yAp = xAp.*-tan(thetaAp);
% (x,y) data for transition with R centered at (xTrOrigin, yTrOrigin)
xTrOrigin = xTrStart + RadiusTr*sin(thetaAp);
yTrOrigin = yTrStart + RadiusTr*cos(thetaAp);
xTr = linspace(xTrStart,xToStart,NUMPOINTS);
yTr = yTrOrigin - sqrt(RadiusTr*RadiusTr-(xTr-xTrOrigin).*(xTr-xTrOrigin));
% (x,y) data for takeoff segment
xTo = linspace(xToStart,xToEnd,NUMPOINTS);
ext = 10; % adding an extension to help takeoffevent trigger
xTo = [xTo + ext];
yTo = yToStart + (xTo-xToStart).*tan(thetaTo);
% stitching together all (x,y) data from inrun
x = [xAp(1:(end-1)), xTr(1:(end-1)), xTo];
y = [yAp(1:(end-1)), yTr(1:(end-1)), yTo];
% evenly spacing stitched (x,y) data
xInrun = linspace(x(1),x(end),500);
xInrunBeg = xInrun(1);
xInrunEnd = xInrun(end);
yInrun = interp1(x,y,xInrun);
% % extended data to use for event
% xExtIn = [xAp(1:(end-1)), xTr(1:(end-1)), xExtTo];
% yExtIn = [yAp(1:(end-1)), yTr(1:(end-1)), yExtTo];
% % evenly spaces data
% xExtInrun = linspace(0,xExtIn(end),500)
% yExtInrun = interp1(xExtIn,yExtIn,xExtInrun);
%%Using ode45 to integrate equations of motions to find velocity of rider along inrun
% reference files: inrun, takeoff_event
dydxInrun = [0 diff(yInrun)./diff(xInrun)];
ddydxInrun = [0 diff(dydxInrun)];
kurvature = ddydxInrun./(1+dydxInrun.^2).^1.5;
tF=20;% final time, seconds
dt=0.01;
tspan=0:dt:tF;
options = odeset('Events',@takeoff_event, 'RelTol', 1e-6);
%states are [x v ]
stateInr_0 = [xRiderStart 1.2]; %initial condition
%integrate ode's 10
%ode=@(t,state) inrun2(t,state);
[ti,stateInr_i,tei,yei,iei]=ode15s(@inrun,tspan,stateInr_0,options);
tstateei=[tei yei iei]
xInr=stateInr_i(:,1);
yInr = interp1(xInrun, yInrun, xInr);
vInr=stateInr_i(:,2);
-----------------------------------------
function [value,isterminal,direction] = takeoff_event(t,state)
global xInrun ext xToStart
x=state(1);
% entry to takeoff ramp
value(1)=x-xToStart;
isterminal(1)=0;
direction(1)=1;
%pause
% extension subtracts from xInrun to actual xTo.
% extension is needed so that event can trigger properly.
% exit from takeoff ramp
xTo=xInrun(end)-ext; % xTo = x takeoff (m)
value(2)=x-xTo;
isterminal(2)=1;
direction(2)=1;
return
-----------------------------------------
function [ stateIRdot ] = inrun(t,stateIR)
% ODEs for inrun acceleration
% state vector = [ x v ]
global xInrun yInrun dydxInrun kurvature
global m CdA mu eta g
x=stateIR(1); % x = horizantal coordinate
v=stateIR(2);% v = velocity along path
dydx=interp1(xInrun,dydxInrun,x);% finding dydx for any value of x
k=interp1(xInrun,kurvature,x);% finding curvature for any value of x
theta=atan(dydx);% angle of slope at any value of x
stateIRdot=[0 0]';% setting initial value of statedot
stateIRdot(1)=v*cos(theta);
N=m*(g*cos(theta)+k*v*v);% normal force
stateIRdot(2)=-g*sin(theta)-eta*v*v-mu*N*sign(v)/m; % since theta<0 sign cancels
txvstatedot=[t x v stateIRdot'];
%pause
end

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

回答（1 个）

Jan
Jan 2013-6-18
To let the integrator determine the event time, the observed variable needs to cross the zero value. To detect this NaN's are not allowed in the evaluated trajectory. So you need a method to define valid values.
Debugging your code would be much easier, if you omit the evil "clear all" call, which deletes all breakpoints also. Compared with this brute impeding of the debugger, it deletes all loaded files from the memory and reloading them wastes a lot of time.

类别

在 Help Center 和 File Exchange 中查找有关 Code Execution 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by