Split Impulses and Joints

[raigan2]

Erin: Awesome! I can't see any position error at all in the demo.

You're right, Baumgarte works great for regular joints; my problem is that I'm trying to get "rigid" joints working (spherical joint constraint + orientation constraint) for use with some physics-driven animation stuff, and Baumgarte was much too bouncy/loose -- even when solving the position+orientation constraints simultaneously (Jacobian is 3x3).

Is the Jacobian always invertible? I see you're using "Damped Newton" (adding a scaled identity to the Jacobian to make sure it's invertible), however your damping parameter ("softness") seems to be set to 0 in the constructor.


Jan: Thanks! I'll compare that with my code and see what the problem is.

[Dirk Gregorius]

@Erin:

You don't consider the contacts in the "super-duper" position correction, but the chain demo looks quite convincing. At least I did not spot any artifacts. Do you see any problems with this?

Obviously this is much more expensive than Baumgarte or split impulses

Isn't this relative? The only other option I see to get more rigid constraints would be to use a direct solver/subsolver. So I wonder if this is a nice option to add this as "post-stabilization" for subsystems that need higher accuracy for position constraints? You also assume small violations of the position constraint. What would be the effect when some position constraints have a bigger error?

Can I extend the mechanism simply to unilateral constraints (e.g. limits)? With contacts I see the problem that we have to call the collision system at least for each outer iteration which then becomes quite expensive I guess...

However, they don't work well with spherical joints.

Can you say why?


Cheers,
-Dirk

PS:
I think there was a typo in the version I just downloaded. In World.cpp line 157 the index should be "j" instead of i...

[Dirk Gregorius]

I made some experiments and switched to your usual solver (SI with Baumgarte and warmstarting). The chain is jitterin quite a lot in my opinion when you use a Baumgarte factor of ~0.2. When I reduce the Baumgarte term to 0.05 the simulation looks much better (more calm). There seems to be a coupling between the number of constraints and the Baumgarte factor. I see that you are playing with CFM (softness) in the suspension demo (Demo . Do you think that we need to adjust these values as a function of the number of constraints in the local system - that is make everything softer for larger system?

Cheers,
-Dirk

[Erin Catto]

raigan: Yes, the Jacobian for the spherical joint is always invertible. I was using softness to simulate a spring damper system for the suspension bridge.

It sounds like you are trying to use motors to play animations. From my experience animations often have large accelerations, requiring massive motor forces. You should expect some smoothing of the motion. However, we should be able to prevent large joint separations.

Dirk: Yes, you can mix and match with position stabilization. You can probably even use post projection along with Baumgarte on the same constraint.

Post-projection doesn't require small position errors to get _a_ solution. Given enough iterations, it should eventually converge to a local solution.

You should be able to use post-projection with unilateral constraints. You will have to accumulate impulses. I don't think warm starting should be used. You could update contact point position and separation without doing collision detection. You couldn't get any new contact points this way, but I think that should not be a problem.

I haven't come to any conclusions on the suspension bridge.

[Dirk Gregorius]

Actually I wonder why we have to iterate at all. Since you know this is a chain (and this works for trees as well) you simply move from the root down to the childs and set the current parent immovable and adjust the child accordingly.

For your particular chain example you simply walk down from top to bottom, set the parent to immovable and correct the position error through only moving the child body. From my experience this damps quite a bit, but I will try this out tomorrow and compare it to the iterative Newton solution. We have some information about the system so we should use it, right?

I also noticed that with any approach this system doesn't come to rest. For Baumgarte this is more noticeable then with post-stabilization. Do you see any reasons why the system doesn't converge against an equlibrium state and comes to rest - that is simply hanging down? Even the small jittering is quite unpretty.

Cheers,
-Dirk

[raigan2]

raigan: Yes, the Jacobian for the spherical joint is always invertible. I was using softness to simulate a spring damper system for the suspension bridge.

It sounds like you are trying to use motors to play animations. From my experience animations often have large accelerations, requiring massive motor forces. You should expect some smoothing of the motion. However, we should be able to prevent large joint separations.

That's odd -- I've been reading a lot about nonlinear solvers (trying to understand the "Position Based Dynamics" paper), and it seems like CFM is exactly the same as some Newton-like methods -- adding a scaled identity matrix to the Jacobian before inversion to make it more well-behaved. I think the ODE docs even suggest that adding some CFM has the same effect -- adding a bit will help probablematic systems converge/avoid blowing up.


Also, AFAIK Newton iteration isn't guaranteed to converge regardless of iteration number -- you have to start close enough to the solution or it may fail to find a minimum. This may be what's happening in my tests; I'm interacting by simply displacing bodies, and this can introduce relatively large position errors.


Finally, I am indeed using the orn part of the "rigid" constraint as a motor, however my problem isn't that the resulting movement is somewhat soft, it's that the joints themselves don't want to converge/a lot of position error is introduced.

For instance, a chain of 3 bodies joined by "rigid" joints should be fixed wrt each other, however in my cases all of the solvers I've tried have failed to keep the joints together. It would be fine if the joint angle isn't met exactly (this is definitely to be expected, I'm just trying to drive the system _towards_ that desired angle), it's that not only are the angles not being met, but in the process the joints themselves are being forced apart -- large position error is present.

If something lands on my gigantic robotic arm, it's fine if the arm isn't able to maintain the same pose -- but it should stay together/bend gracefully, instead of large gaps appearing at the joints! Currently I'm just trying to get an arm suspended in space to flex/unflex nicely, but the joints are either coming apart or completely exploding as soon as something (even gravity) tries to oppose the flexion.

[Erin Catto]

Dirk: right, but I'm presenting a general solver. As soon as we get a loop, or include contact in the post-projection, we will need to iterate. However, I'm not against exploiting special cases.

Semi-implicit Euler is symplectic, so it doesn't damp numerically. The damping we get is from the constraint solver, but we shouldn't depend on that to add numerical damping in all cases. I added some angular damping and the chain fell asleep. If I was putting a chain like this in a game and wanted it to sleep quickly, I would put friction motors on each joint (target velocity zero).

[Erin Catto]

raigan: You are correct about CFM. The Hessian diagonal is often increased in optimization in order to drive the solution towards a local minimum. Otherwise the solution can get stuck at a local maximum or saddle point.

I think that people used to use CFM in ODE with the old direct solver. Direct solvers are notorious for numerical sensitivity and the diagonal adjustment using CFM helps to improve the numerics. One thing that causes direct solvers a lot of grief are systems with redundant constraints because the Jacobian has linearly dependent rows.

The ODE quickstep solver uses Projected Gauss-Seidel (PGS) which is roughly equivalent to Sequential Impulses (SI). Both PGS and SI are free of numerical sensitivity. They handle redundant constraints without any problems. Of course the convergence of PGS and SI is slow.

CFM also has another use. That is to simulate compliant joints. So you can use it to model stretchy ropes or vehicle suspensions. I posted an update to Box2D that correctly uses CFM (softness) on the suspension bridge, Demo 8:

http://www.gphysics.com/files/Box2D_SuperSplit.zip

I've done a lot of work with Newton solvers for multi-body linkage systems. They have always converged for me, but you can get into situations where you get the wrong solution (the linkage flips) or you get stuck at a local maximum or saddle point. So, in my experience, Newton solvers can handle quite a bit of error.

I'm surprised at your difficulties. I think your goal should be achievable with existing methods. If you'd like to give me a simple test case that fails in your setup, I can try to emulate it in Box2D.

[Jan Bender]

Can anybody tell me how fast the nonlinear solvers like Newton iteration converge. My method uses an approximation to get a linear equation which can be easily solved. Due to this approximation you don't get the exact solution in the first iteration but usually after two iterations (I made several tests with this) the error is smaller than 10^-6 m. So I am interested if the nonlinear solvers can beat this or not.

Jan

[Dirk Gregorius]

I don't see a big difference between your approximation and the Newton solver. Actually I would argue that you solve the linear system in the inner loop using an Gauss-Seidel like iterative approach with only one outer iteration. The only difference is that you adjust the velocity such that position constraints will be satisfied at the end of time step, but principally I would argue that this is the same.

A reason for the fast convergence (note that this is only a guess) might be that you also correct the velocity constraints so you have a pretty good approximation in the first place. What is questionable in your approach is if it is physical correct to relate position constraints and velocities? As Erin pointed out the correct solution would be:

1) Solve for the Langrange multipliers (solving a linear system)
2) Integrate the velocity using whatever method you want
3) Project the integrated velocities onto the (velocity) constraint manifold solving a linear system
4) Integrate the positions
5) Project the integrated positions onto the (position) constraint manifold solving a non-linear system (using e.g. a Newton method since you can expect that the solution is close to the correct global solution).

Cheers,
-Dirk

[Dirk Gregorius]

@Erin:

I have two questions:

a) Can you explain/describe a bit what you did for the suspension bridge, e.g. how you choose the stiffness and damping. Maybe you can also explain why this works and what made the bridge bouncy without CFM?

b) I played also a little with the chain example (Demo0). Using the "super-duper" stabilization gives me a very nice chain. When I set the splitImpulse variable in the joint to "false" the chain gets very bouncy/jittering again. Could the behaviour of the chain be improved using simple Baumgarte stabilization (maybe using some softness/CFM? What is the reason for the jittering - does the stabilization term add to much memory back to the system and can this be avoided?


Great work!
-Dirk

[Jan Bender]

A reason for the fast convergence (note that this is only a guess) might be that you also correct the velocity constraints so you have a pretty good approximation in the first place. What is questionable in your approach is if it is physical correct to relate position constraints and velocities?

But I don't correct the velocity constraints in the same iteration loop. It converges also fast if I don't correct the velocities at all.

I have written a whole paper with my professor concerning the physical correctness of the method. This paper also contains a proof for the physical correctness. So I wouldn't say that it is questionable.

Jan

[raigan2]

I've done a lot of work with Newton solvers for multi-body linkage systems. They have always converged for me, but you can get into situations where you get the wrong solution (the linkage flips) or you get stuck at a local maximum or saddle point. So, in my experience, Newton solvers can handle quite a bit of error.

I've been trying to get a particle-based angle constraint working (following "Position Based Dynamics") and while it often converges, if you take one end of a chain and yank it, things tend to explode.

I'm surprised at your difficulties. I think your goal should be achievable with existing methods. If you'd like to give me a simple test case that fails in your setup, I can try to emulate it in Box2D.

I've just quickly implemented orientation constraints into your latest demo, it's uploaded here: http://www.harveycartel.org/raigan/Box2D_r.zip

Part of my problem I think, is that I've been using very low numbers of iterations, but even when using 10 iterations with supersplit (5 inner iterations as in your code), the same problems happen.

Demo 1 is just some boxes rigidly joined to each other.

Demos 2/3/4, 5/6/7, and 8/9/0 are 3-,4-, and 5-link chains respectively, with 22.5deg orientation constraints between each link.

The three versions of each chain show [baumgarte joints + baumgarte orn], [supersplit joints + baumgarte orn], and [supersplit joints + supersplit orn] respectively.


Supersplit is definitely an improvement _when_ it works (demo 4), however there are cases where it performs much worse (demo 0). Baumgarte is unfortunately much too "jiggly" for use with animation -- it's fine if things converge slowly/are soft, but springy oscillations don't work.

My orn constraint is seperate at the moment, I have a joint+orn constraint which I'm also going to try adding to box2D, however I suspect I'll have the same problems (seen in demo 0).

My main concern is keeping the joints together; I don't care if the orientations aren't met exactly, or are met slowly, but they can't tear the joints apart as they're currently doing!

Anyway, thanks for the input and help, I don't really understand what's going wrong in demo 0.

oh, also -- I had to disable warm-starting for the orn constraints or things would blow up immediately.

[Erin Catto]

Dirk: I will write up something on the bridge when I have time.

It is hard to say exactly why Baumgarte can sometimes lead to instabilities. However, it is easy to test if Baumgarte is to blame. Just set the bias factor to zero and run with no position correction. You might be surprised how well things hang together.

We should not be upset if Baumgarte leads to instabilities. It is, after all, a hack. Post-projection is the most robust method.

raigan: It seems like you ran straight towards one of the hardest problems for iterative solvers. You are trying to create rigid structures using weld joints. Creating cantilever beams out of multiple bodies and weld joints is quite challenging, especially under the load of gravity.

Iterative methods don't do well with large unbalanced torques. For these cases I usually create counter balance ghost bodies (with no collision). If your application is not for serious engineering, then tweaking the model is often the best course of action.

If my goal was to create a physics engine capable of simulating long motorized chains under gravity, then I would use a joint coordinate model (often called Featherstone).

I played with your demos and eventually wrote a fixed joint that seeks a target rotation. It solves all three joint constraints at once using a 3x3 matrix. It works okay with post-projection, but there are some oscillations. Baumgarte blows up, except for very small bias factors. You can check it out here:
http://www.gphysics.com/files/Box2D_FixedJoint.zip

[Dirk Gregorius]

If my goal was to create a physics engine capable of simulating long motorized chains under gravity, then I would use a joint coordinate model (often called Featherstone).

For my knowledge Featherstone doesn't work woth motors. Do you have some reference for this?