ODE OutputFcn returning fewer points than expected

3 visualizzazioni (ultimi 30 giorni)
Hi,
I'm using OutputFcn to return exta information from ode45 and whilst the returned signal looks as I expected it to it's only comprised of 93 points rather than 1001.
tvec is my tspan vector and contains 1001 points for a run of 0 to 5 seconds at 200 Hz.
The ODE setup is:
% Set up ODE function to pass to solver
odeFcn = @(t, x) Three_DOF_Half_Car_ODE(t, x, ix, iz, ip, tvec, mb, Iby, lf, lr, h, ...
kfz, cfz, kfx, cfx, adf, adf_dbz, krz, crz, krx, crx, adr, adr_dbz, Fxr_max);
% Set up options for output function
options = odeset('OutputFcn',@(t, x, flag) myOutputFcn(t, x, flag, ix, tvec, lr, ...
krx, crx, adr, adr_dbz, Fxr_max));
% Run solver
[T,x] = ode45(odeFcn, tvec, x0, options);
The intention is that Fxr is calculated during integration and its value capped under certain conditions (I'm aware this is an abuse of the ODE solver). I want Fxr out so I repeat the same equations in the output function, which contains:
function status = myOutputFcn(t, x, flag, ix, tvec, lr, krx, crx, adr, adr_dbz, Fxr_max)
persistent Fxr_out
ix = interp1(tvec, ix, t); % interpolates ix at t within tvec
switch flag
case 'init'
tanadr = tan(deg2rad(adr_dbz.*(x(1)+lr.*x(3))+adr)); % anti-dive tangent wrt body displacement
xir = (x(1)+lr.*x(3)).*tanadr; % contact patch longitudinal displacement...
xdotir = (x(4)+lr.*x(6)).*tanadr; % ..and velocity
Fxr = krx.*(-xir-x(2))+crx.*(-xdotir-x(5)); % rear axle longitudinal force
if ix ~= 0 % if longitudinal acceleration is negative...
if Fxr < Fxr_max % if Fxr exceeds limit...
Fxr = Fxr_max; % cap Fxr to Fxr_max
end
end
Fxr_out = Fxr;
case ''
tanadr = tan(deg2rad(adr_dbz.*(x(1)+lr.*x(3))+adr)); % anti-dive tangent wrt body displacement
xir = (x(1)+lr.*x(3)).*tanadr; % contact patch longitudinal displacement...
xdotir = (x(4)+lr.*x(6)).*tanadr; % ..and velocity
Fxr = krx.*(-xir-x(2))+crx.*(-xdotir-x(5)); % rear axle longitudinal force
if ix ~= 0 % if longitudinal acceleration is negative...
if Fxr < Fxr_max % if Fxr exceeds limit...
Fxr = Fxr_max; % cap Fxr to Fxr_max
end
end
Fxr_out = [Fxr_out, Fxr];
case 'done' % when it's done
assignin('base','Fxr',Fxr_out); % get the data to the workspace.
end
status = 0;
As i say, the resulting data in Fxr is correct. However, it only contains a fraction of the expected number of points. Why is this?
Regards,
Simon.

Risposte (1)

Jan
Jan il 26 Gen 2021
Modificato: Jan il 26 Gen 2021
The code looks okay. Check again if tspan has really 1001 time points.
Use the debugger to check, what's going on. Set a breakpoint in the output function and see, when it is triggered.
I'd stay away from using assignin to create a variable in the base workspace. What about:
function status = myOutputFcn(t, x, flag, ix, tvec, lr, krx, crx, adr, adr_dbz, Fxr_max)
persistent Fxr_out
...
status = 0;
switch flag
case 'init'
...
case ''
...
case 'done'
% Nothing
case 'flush'
status = Fxr_out;
Fxr_out = [];
end
end
Then calling this function after the integration to flush the internal persistent buffer is save, because the data are actively requested. Poking it by assignin into the base workspace is susceptible as all global access to variables: If there is another forgotton assignin anywhere in your code, the confusion is perfect.
Note: tand() is nicer than tand(deg2rad()).
  4 Commenti
Simon Aldworth
Simon Aldworth il 29 Gen 2021
Hi Jan,
I have managed to reformat my code to deal with multiple values in t:
switch flag
case 'init'
tanadr = tand(adr_dbz.*(x(1,:)+lr.*x(3,:))+adr); % anti-dive tangent wrt body displacement
etc...
I've also changed the if statement to deal with multiple values in ix and Fxr:
idx = Fxr < Fxr_max & ix~=0;
Fxr(idx) = Fxr_max;
So, thanks for your help and guidance.
Just one last question: I can't find any help documentation for using 'flush' with flag. Where did that come from and is there a help page hidden somewhere?
Thanks and regards,
Simon.
Simon Aldworth
Simon Aldworth il 5 Feb 2021
@Jan Hi, can you tell me about using 'flush' with flag - is there any documentation for that?
Thanks.

Accedi per commentare.

Categorie

Scopri di più su Programming in Help Center e File Exchange

Tag

Prodotti


Release

R2018b

Community Treasure Hunt

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

Start Hunting!

Translated by