Physical realism of damped pendulum

[ratatouille]

Hey all,

I am comparing the physical accuracy of a simple pendulum with a damper and torsional spring. My benchmark is an implementation of the same dynamical system in Python. The problem is that the results seem to diverge quite a bit.



I am using the btGeneric6DofSpring2Constraint. I have already played with CFM/ERP (on all 6 axes) which does not seem to many any difference. Linear and angular damping on both objects (the pendulum cube and the cube at the point of rotation) are 0. maxSubSteps is equal to 0 and the fixed simulation timestep is 5 ms. Gravity is 0.

I know that Bullet is not intended to be super physically accurate, but I'm sure we can get closer than this. Where else could I look?

[drleviathan]

Does your python implementation compute the theoretical curve or does it do a simulation with timesteps?

If it computes the theoretical curve then one thing to try is to play with the timestep of the bullet simulation to see if it get better with smaller steps and worse with larger.

If it performs a simulation with timesteps then both simulations are suspect and you need to plot the theoretical curve to see how each compares.

[Basroil]

Are you using MLCP solver or SI? If SI, are your iteration steps high enough?

A swinging pendulum is non-linear, so you should expect some error using linear or partially-linearized systems. How much error is involved depends on a lot of factors (though I'm sure you already knew that)

[ratatouille]

I tried a bunch of different solvers: Lemke, Dantzig, Gauss-Seidel and the SI solver with a constant timestep and no substepping. None make any visual difference to the plot. I also reduced the timestep to 0.1 ms, again no visual difference. That's why I am thinking it is not the solver, but some hidden damping parameter or something that I'm missing.

Some further things ruled out: additionalDamping is false, activationState is DISABLE_DEACTIVATION.

thanks for any leads...

[Basroil]

I tried a bunch of different solvers: Lemke, Dantzig, Gauss-Seidel and the SI solver with a constant timestep and no substepping. None make any visual difference to the plot. I also reduced the timestep to 0.1 ms, again no visual difference. That's why I am thinking it is not the solver, but some hidden damping parameter or something that I'm missing.

Your "reduced the timestep" means what? Did you change the time between calls to stepSimulation (first parameter), or change substep and simulation framerate (second+third parameters)? If all you changed is the first, then of course the solution won't change!

[ratatouille]

My call looks like this:

Code: Select all

int maxSimSubSteps = 0;		// constant timestep
int numSimSteps = m_dynamicsWorld->stepSimulation(dt, maxSimSubSteps);

It is my understanding that for maxSimSubSteps = 0, the step of size dt is not further subdivided. (See http://bulletphysics.org/Bullet/BulletF ... d8f8dd6294 )

[Basroil]

My call looks like this:

Code: Select all

int maxSimSubSteps = 0;		// constant timestep
int numSimSteps = m_dynamicsWorld->stepSimulation(dt, maxSimSubSteps);

It is my understanding that for maxSimSubSteps = 0, the step of size dt is not further subdivided. (See http://bulletphysics.org/Bullet/BulletF ... d8f8dd6294 )

I suggest you take a look at http://bulletphysics.org/mediawiki-1.5. ... _the_World instead. Your code is not changing the resolution properly, so you need to add the third term.

[ratatouille]

I've also tried a call of the form:

Code: Select all

int maxSimSubSteps = 1;		// constant timestep
int numSimSteps = m_dynamicsWorld->stepSimulation(dt, maxSimSubSteps, dt);

with no changes (dt = 1 ms or 0.1 ms). I've also upgraded the forward Euler solver in the Python validation script to a scipy explicit Runge-Kutta ('dop853'), with again no changes.

Looks like additional damping term still hidden in Bullet somewhere.

[drleviathan]

Yes, when maxSimSubSteps = 0 Bullet is supposed to accept the provided dt at face value and only simulate forward the specified amount of time. The reason this is not recommended for real-time games is because if you feed it fixed timesteps then your simulation will not run at real-time, and if you feed it variable timesteps then you risk instability if you happen to give it too large of a timestep, and will get somewhat irreproducible results even if you happen to always give it small timesteps.

For your application the simulation doesn't have to run at real-time so supplying small fixed timesteps is the correct way to go.

Assuming you were providing dt in units of seconds rather than milliseconds, and were getting the same results for 1/60th sec and 1/10000th sec, and across all of the bullet solvers, then either: (a) there really is some damping in bullet that is independent of the solver or (b) there is a problem with the yardstick you're comparing them to.

[ratatouille]

OK, it looks like there's some funny stuff going on with the btGeneric6DofSpring2Constraint. Plotting the spring torque (via setJointFeedback->m_appliedTorqueBodyA) reveals discontinuities every time after crossing the resting position:



The discontinuity is not a single jump, but exists for several timepoints (each timepoint = 1 marker circle):

[drleviathan]

I find this unsurprising. For stability and optimization we often add little thresholds that disable behavior when velocities are near zero. Often it is one of those tradeoffs between computation speed and correctness.

The lesson here might be: Bullet is probably not the best physics engine for performing very accurate simulations of physics interactions, but is good for providing real-time simulations that pass an eyeball test (I doubt that glitch would be visible to someone watching it in real-time).

[Basroil]

There's a value called D6_LIMIT_ERROR_THRESHOLD_FOR_ROTATION, checked to see if it changes your issue?

Torque goes to zero at the minimum/maximum of a pendulum (assuming 360 degree movement, I guess just minimum in your case?), where you have two issues I can think of:
1) your torque becomes small enough to get lost in a 32bit float (not as likely I guess)
2) the pendulum accumulates enough joint error to cause heavy handed corrections, since torque=0 also means a=g+a bit (rather than g-a bit) if we are forced to use old estimates for a

Do you see a similar jump in position/velocity?

[ratatouille]

Some good news and bad.

The jump in spring torque has been eliminated by changing lines 824/825 in btGeneric6DofSpring2Constraint.cpp:

Code: Select all

-                       info->m_lowerLimit[srow] = -maxf > 0 ? 0 : -maxf;
-                       info->m_upperLimit[srow] = -minf < 0 ? 0 : -minf;
+                       info->m_lowerLimit[srow] = -maxf > 0 ? 0 : -f;
+                       info->m_upperLimit[srow] = -minf < 0 ? 0 : -f;

Bad news is that the change in dynamics is virtually nil. Bullet vs. Python (Runge-Kutta) solution curves still look like image in 1st post.

[ratatouille]

Sorry, forgot to mention. Basroil, thanks for your suggestions. D6_LIMIT_ERROR_THRESHOLD_FOR_ROTATION makes no apparent difference. I also recompiled everything with double precision floats. Again no evident change.

[ratatouille]

Could angular inertia be the culprit? In the Python case, the pendulum is ideal, with a point mass. In Bullet, the mass is distributed over a rigid body of nonzero size (i.e. the swinging body of the pendulum). For a long enough pendulum and a low enough mass (L = 10 m, m = .1 kg) the results overlap.

I tried setting the inertia tensor and/or the angular factor to 0 to try and disable this effect. But when I do this, the spring constraint doesn't work any more and the pendulum just flies off to the right.