One column in return from ode45 is all zero
조회 수: 1(최근 30일)
표시 이전 댓글
Simon Aldworth 2021년 2월 11일
댓글: Star Strider 2021년 2월 12일
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
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))-...
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))-...
dx(7) = x(8);
dx = dx.';
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.
댓글 수: 0
Star Strider 2021년 2월 11일
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.
댓글 수: 8
Star Strider 2021년 2월 12일
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.
Find more on Programming in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!Start Hunting!