One column in return from ode45 is all zero

Hi all,

I've setup a simple 3 DOF model of a car to examine its pitch, bounce and longitudinal behaviour during braking events. To limit the rear brakes in a somewhat realistic way I have included a cap on the rear contact patch longitudinal force and use the event location function called from ode45 to stop and restart the integration each time the cap is crossed. To ensure that the jacking forces from the rear suspension are effective, I've added a fourth DOF to capture the motion of the grounded end of the spring-damper that produces the rear longitudinal forces. When the rear braking forces hit the allowable limit and become capped the grounded end of this spring-damper can move so that when the cap is removed the new position for this grounded part will cause the vehicle to remain in a somewhat pitched (jacked) position. This is perhaps easier understood in the following image where all springs also include dampers and the rear horizontal spring-damper is connected to a friction element that will allow the grounded end of the spring to move once a certain load is achieved.

The model returns seemingly plausible results with the exception of the velocity of the grounded part. This is zero throughout, despite the fact that all the other signals returned are non-zero.

My ode equations are as follows:

function dx = 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, EBD_state)

% equations of motion
% x(1) = body vertical displacement
% x(2) = body longitudinal displacement
% x(3) = body pitch displacement
% dx(1) = x(4) = body vertical velocity
% dx(2) = x(5) = body longitudinal velocity
% dx(3) = x(6) = body pitch velocity
% dx(4) = body vertical acceleration
% dx(5) = body longitudinal acceleration
% dx(6) = body pitch acceleration
% x(7) = rear contact patch ground displacement
% dx(7) = x(8) = rear contact patch ground velocity
% dx(8) = rear contact patch ground acceleration

ix = interp1(tvec, ix, t); % longitudinal acceleration
iz = interp1(tvec, iz, t); % not in use
ip = interp1(tvec, ip, t); % not in use

tanadf = tand(adf_dbz.*(x(1)-lf.*x(3))+adf); % front anti-dive angle wrt body
tanadr = tand(adr_dbz.*(x(1)+lr.*x(3))+adr); % rear anti-dive angle wrt body

xif = -(x(1)-lf.*x(3)).*tanadf; % front contact patch displacement 
xir = (x(1)+lr.*x(3)).*tanadr; % rear contact patch displacement 
xdotif = -(x(4)-lf.*x(6)).*tanadf; % front contact patch velocity 
xdotir = (x(4)+lr.*x(6)).*tanadr; % rear contact patch velocity 

% Set contact patch ground velocity and acceleration to initial value of zero
% x(8) = 0; % NOT SURE IF THIS IS NECESSARY
dx(8) = 0; % rear contact patch ground acceleration is zero unless modified below

Fxr = krx.*(-xir-x(2)+x(7))+crx.*(-xdotir-x(5)+x(8)); % Rear axle longitudinal force

% Switch value of Fxr depending on EBD_state
if EBD_state == 1 % rear brake force limit is active
    if Fxr < Fxr_max % if rear brake forces exceeds limit (negatively)
        Fxr = Fxr_max; % limit rear force to maximum (negative) value
        x(8) = xdotir; % set rear contact patch ground part velocity to match contact patch
        dx(8) = (dx(4)+lr.*dx(6)).*tanadr; % set rear contact patch ground acceleration to match contact patch
    end
end

dx(1) = x(4);
dx(2) = x(5);
dx(3) = x(6);
dx(4) = 1./mb.*(kfz.*(-x(1)+lf.*x(3))+cfz.*(-x(4)+lf.*x(6))+krz.*(-x(1)-lr.*x(3))+crz.*(-x(4)-lr.*x(6))-...
    (kfx.*(-xif-x(2))+cfx.*(-xdotif-x(5))).*tanadf+Fxr.*tanadr+mb.*iz);
dx(5) = 1./mb.*(kfx.*(-xif-x(2))+cfx.*(-xdotif-x(5))+Fxr+mb.*ix);
dx(6) = 1./Iby.*(-kfz.*lf.*(-x(1)+lf.*x(3))-cfz.*lf.*(-x(4)+lf.*x(6))+krz.*lr.*(-x(1)-lr.*x(3))+crz.*lr.*(-x(4)-lr.*x(6))-...
    kfx.*h.*(-xif-x(2))-cfx.*h.*(-xdotif-x(5))-Fxr.*h-...
    lf.*(kfx.*(-xif-x(2))+cfx.*(-xdotif-x(5))).*tanadf+lr.*Fxr.*tanadr+Iby.*ip);
dx(7) = x(8);

dx = dx.';
end

I know from setting breakpoints and from the returned results that the above 'if' statement is used and that during the solve x(8) is given a non-zero value from xdotir. However, the results given back to my main model code have all elements of x(8) equal to zero.

I have avoided putting in too much code for now to keep this relatively short and to get a check to see if the above code makes sense. However, I'm happy to add more info where needed.

Thanks,

Simon.

Answer 1

[Star Strider]

Assuming that ‘x(8)’ is 0 and not simply a very small value, There’s obviously something about the nested if blocks that’s not working.

First, check to see whether ‘EBD_state’ is exactly equal to 1. If it isn’t, no matter how close it is to 1, that equality will never be true.

Second, remove the trailing semicolon from any of the assignments inside the inner if block. I suspect that it’s not being entered, so ‘dx(8)’ never changes. Then, figure out what the problem is.

Also, note that the MATLAB (and likely all) numeric ODE integration routines do not handle discontinuities well, so that if ‘dx(8)’ has any sort of ‘significant’ discontinuity (not differentiable), the result may not be reliable.

[Simon Aldworth]

Hi, thanks for the quick response.

Taking your points in order:

EBD_state is exactly one; I set it so using

EBD_state = 1;

Inserting EBD_state as a line inside the 'if' statement return EBD_state = 1 in the output window.

Removing the trailing semi-colon from the line setting x(8) to equal xdotir returns multiple instances of x in which x(8) is non-zero, e.g.:

x = 8×1    
   -0.0225
    0.1061
    0.0390
   -0.0724
   -0.2188
    0.2208
    0.0235
    0.1759

However, removing the semi-colon from the dx(8) assignment returns the following (interlaced with EBD_state and x(8) returns):

EBD_state = 1
x = 8×1    
   -0.0112
    0.0873
    0.0099
   -0.1019
    0.4745
    0.1714
    0.0024
    0.0959
dx = 
     0     0     0     0     0     0     0     0

Could the problem be the different arrangement of the two vectors x and dx?

I asked for the same information from the first run of the solver (by inserting EBD_state, x and dx just before the 'if' statement), where EBD_state is zero and the 'if' statement ignored and I got the same results:

EBD_state = 0
x = 
   -0.0032
    0.0372
    0.0011
   -0.0941
    0.5530
    0.0334
         0
         0
dx = 
     0     0     0     0     0     0     0     0

It seems that at this point dx is always zero. If this is the case, how should I formulate my statement for dx? Should I use an alternative like:

dx(8) = ((gradient(x(4))/gradient(t))+lr.*((gradient(x(6))/gradient(t))).*tanadr

Regards, Simon.

[Star Strider]

My pleasure!

I’m currently willing to believe that the values may not actualkly be 0, simply very small. To rule that out, change to:

format long

or:

format longE

and then see if they are all still 0.

Also, ‘x(8)’ (and the others) are provided to ‘Three_DOF_Half_Car_ODE’ by the ODE integration function you’re using. Don’t change it in your code.

[Simon Aldworth]

Good morning. Here are the results after formatting to longE:

EBD_state = 
     1
x = 
    -1.732958657425211e-02
     1.078626330179621e-01
     2.354065970906365e-02
    -8.688544526111880e-02
     1.516750230399420e-01
     2.392392784781251e-01
     1.156031091187122e-02
     1.757679521402903e-01
dx = 
     0     0     0     0     0     0     0     0

It seems that interogating dx at this point in the code is fruitless, but since the the model runs and the parts certainly do accelerate, then my main concern is why x(8) is non-zero here and yet zero when returned to the model.

Regards, Simon.

[Star Strider]

Good morning!

Note that I continue to maintain that re-defining ‘x(8)’ in your code is a terrible idea. It should be whatever the ODE integrating function defines it as, and sends it to ‘Three_DOF_Half_Car_ODE’. This is going to cause problems and will not produce the result you want. Calculating with the arguments is appropriate. Re-defining them is not.

That aside, since the if block logic seems to be working correctly, my guess is that in the ‘dx(8)’ assignment, either all the components inside the parentheses or ‘tanadr’ (or all of them) are 0.

[Simon Aldworth]

Hi, I'm so sorry but I think I just took on board the point you're making about re-defining x inside the ODE. I thought you were referring to something else, until it just twigged! My apologies.

So, taking a step back: the content of the 'if' statement is only triggered when there is a re-start of the integration, as the event location feature will have halted the integration when Fxr crosses Fxr_max. I think this means that I can place ALL the dx statements inside the 'if' statement, with one set being used if EBD_state ==1 and Fxr < Fxr_max, and a second set for all other times. Is this better practice?

If not, what is best practice for running two sets of differential equations that are used when certain conditions are met?

A problem I think I might have if I took my suggested route is that I am calculating Fxr inside the ODE and so it relies on knowing the current state of the model. How does one know the current state of the model if the integration for the current moment hasn't taken place?

Regards.

[Star Strider]

NMo worries! Just note that re-definig them could cause problems (and make the solution garbage).

The if logic doesn’t appear to be the problem, since it appears to trigger approipriately.

I have no idea what you’re actually doing, so I’m just helping you troubleshoot the problem you posted. If you believe that putting all the derivaitve calculations inside the if block is appropriate, then do so. Note that switching between the two sets of ‘dx’ calculations could cause problems if the transition between them is not smooth (differentiable as opposed to non-differentiable). If that’s relatively continuous, that shouldn’t be a problem. (This is the sort of situation the ‘Events’ function solves, and that you apparently use correctly for other discontinuities. You may simply need to experiment and see if the result makes sense.)

[Simon Aldworth]

Ok, thanks.

Putting all the equations as a block inside the 'if' statement has improved things - I now get sensible values back for all elements of x. Unfortunately, it's also highlighted a flaw in the model so I'm off to work on that now.

Do you happen to know the answer to the question about calculating Fxr inside the ODE? If Fxr is a function of elements of x, then is the result based on the last known condition of the model or is it wrapped up into the integration at the current value of t?

I have seen similar cases on here where an if statement depends on the value of x and dx is adjusted accordingly.

Thanks for your help; I'll accept your response as the answer.

Regards.

[Star Strider]

As always, my pleasure!

‘How does one know the current state of the model if the integration for the current moment hasn't taken place?’

The current state is the previous integration result, so ‘Fxr’ is being calculated appropriately. The ‘dx’ vector will define the next state.

‘I have seen similar cases on here where an if statement depends on the value of x and dx is adjusted accordingly.’

That depends on the model being solved. It may be appropriate to include the current values of the derivatives in calculations subsequent to their being defined. In any event, the derivatives are defined using the current values of the variables in order to calculate the next state.