ODE OutputFcn returning fewer points than expected

Question

Simon Aldworth 2021-1-26

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/727703-ode-outputfcn-returning-fewer-points-than-expected

评论： Simon Aldworth 2021-2-5

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Jan 2021-1-26

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/727703-ode-outputfcn-returning-fewer-points-than-expected#answer_606913

编辑：Jan 2021-1-26

在 MATLAB Online 中打开

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 个评论
显示 2更早的评论隐藏 2更早的评论

Simon Aldworth 2021-1-27

编辑：Simon Aldworth 2021-1-27

Hi Jan,

Thanks for answering and for the tips.

I set a breakpoint in the output function and was very surprised at the results. I had expected that t would be a single value representing the most recent successful integration time step. However, what I see appears to be otherwise.

The first loop through gives the following variable values and statement returns:

tvec = [0,0.005,0.01,0.015 etc] (1x1001 double)

flag = 'init'

t = [0,5]

ix = [0,0]

Fxr = 0

if ix ~= 0 is FALSE

Fxr_out = 0

It appears to me that t is just showing the tspan range and so ix interpolates accordingly. Is this presumably just a response to the initialisation loop? I don't know how the if statement reacts to ix in vector format. As it happens, Fxr and Fxr_out contain the expected result. tvec is as expected.

Moving on to the second loop, the above list changes to:

tvec = [0,0.005,0.01,0.015 etc] (1x1001 double)

flag = ''

t = 0.005

ix = 0

Fxr = 0

if ix ~= 0 is FALSE

Fxr_out = [0,0]

This is entirely as expected. However, things become problematic on the third loop through:

tvec = [0,0.005,0.01,0.015 etc] (1x1001 double)

flag = ''

t = [0.0100,0.0150,0.0200,0.0250,0.0300]

ix = [0,0,0,0,0]

Fxr = 0

if ix ~= 0 is FALSE

Fxr_out = [0,0,0]

Now t is a vector of 5 values, rather than the expected single value (I assumed it would be just the first one, 0.01). I won't post any further loops, except to show that t then starts from the LAST element value from the previous loop:

t(1) = 0.0350

With t skipping forward like that it is little wonder I only get 93 points back. But why is this happening?

Many thanks,

Simon.

Simon Aldworth 2021-1-29

在 MATLAB Online 中打开

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 2021-2-5

@Jan Hi, can you tell me about using 'flush' with flag - is there any documentation for that?

Thanks.

请先登录，再进行评论。

ODE OutputFcn returning fewer points than expected

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（1 个）

4 个评论
显示 2更早的评论隐藏 2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

ODE OutputFcn returning fewer points than expected

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

回答（1 个）

4 个评论 显示 2更早的评论隐藏 2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论