Position Based Dynamics (PBD)

Please don't post Bullet support questions here, use the above forums instead.
Post Reply
zzzzrrr
Posts: 10
Joined: Mon Apr 16, 2007 7:11 pm
Location: Washington, DC

Position Based Dynamics (PBD)

Post by zzzzrrr »

Hi, I've been reading through Muller's paper, and I'm interested in applying PBD to 2D rigid bodies. Does anyone have a suggestion in regard to updating angular velocity and applying friction / restitution?

Thanks!
Dirk Gregorius
Posts: 861
Joined: Sun Jul 03, 2005 4:06 pm
Location: Kirkland, WA

Re: Position Based Dynamics (PBD)

Post by Dirk Gregorius »

For angular constraints you check the current relative angle and set it back to the initial relative angle. Then as for the linear velocities you set the angular velocity as the difference between the old and the new angle divided by the time step.

I am not aware of any good solutions for friction or restitution in position based dynamics since you don't have any geometric (position) constraints for this. Several cheats exist. You can search the forum for some suggestions. I think Mueller also suggests something for friction in his paper.

-Dirk
raigan2
Posts: 197
Joined: Sat Aug 19, 2006 11:52 pm

Re: Position Based Dynamics (PBD)

Post by raigan2 »

We have a 2D rigid-body simulator based on Jakobsen/Muller, however each body modeled as is a pair of particles connected by a stick rather than a true 2D rigid bodies (i.e 2D position + 1D orientation). You should definitely be able to apply the same ideas (least-squares fitting of the system state to the constraint/error functions via Gauss-Newton) to proper 2D rigid bodies.

Our "stick-based" approach was chosen in order to simplify the calculation of the Jacobians; this is similar to Jakobsen's idea of using a triangle of particles, but simpler and faster. You can of course embed constraint geometry (collision shapes, joint positions or directions/axes) in each body's space.

To support velocity constraints, there are a couple options:

1) Explicitly solve at both position and velocity level. You can just plug in velocities instead of positions into the solver, for the most part the Jacobian for positions is the same as for velocities. For example to constrain the position of two points p0,p1 along an axis parallel to n you would have C = Dot(n,p1 - p0), to constrain the velocity of the points you just use C = Dot(n,v1 - v0). In this case the velocity-level Jacobian is the same as for the position constraints, except that you would replace p with v in the formula.

I'm not sure if it's best to solve positions and then velocities, or vice-versa. We used velocity-before-position because this makes it like semi-implicit Euler integration, and also Box2D uses this order so it must be good :)


2) Implicitly solve velocity constraints at the position level: since you know that the velocity v will be inferred as: v = p_post - p_pre (where p_pre is the old/before-solving position and p_post is after solving), you can infer how to alter p_post in order to satisfy a velocity constraint. For the above example, you would have C = Dot(n, (p1_post - p1_pre) - (p0_post - p0_pre)), where the _pre values are taken as fixed/constant and you solve for the _post values.

Of course, this will violate some position constraints (since you're moving the position of things), but that's no different than any other situation -- through iterative relaxation the solver finds an amicable solution which violates all constraints the least amount possible. The benefit of this approach is that you're only solving at the position level so it's less work. When we first thought of it I expected problems since, for example with friction, you'll be moving the contact point along the surface tangent, which will move it out of collision if the surface isn't flat (i.e if it's two circles colliding for instance). In practice we haven't noticed this to be a problem, however we're not doing really hard-core accurate simulation.


Both approaches should let you model friction.

For constraints related to orientation we define directions as vectors and use atan2, i.e C = atan2(d1) - atan2(d0) would constrain d0 and d1 to be parallel; angular velocity is inferred from pre-solver state as with linear velocities, i.e C = (atan2(d1_post) - atan2(d1_pre)) - (atan2(d0_post) - atan2(d0_pre)) would constrain d1 and d0 to rotate at the same rate.

A more interesting/complicated problem is how to properly model motors in a position-based simulation. We've tried two different methods, one is to just copy what Box2D does (i.e motors are constraints whose summed response is clamped to within some upper/lower limit) and another is much weirder; neither is perfect so far unfortunately.

Collision is also a pain in comparison to velocity-level solvers such as Box2D; since the solver can move positions, you really should run collision detection at each solver iteration, but this is quite expensive.
zzzzrrr
Posts: 10
Joined: Mon Apr 16, 2007 7:11 pm
Location: Washington, DC

Re: Position Based Dynamics (PBD)

Post by zzzzrrr »

Many thanks for the great insight! Using a stick model is a very interesting concept, although I'm attempting to take the traditional RB approach before I try anything fancy :D
raigan2 wrote: I'm not sure if it's best to solve positions and then velocities, or vice-versa.
Since the solver is updating position, why also constrain linear velocity? Why not update linear velocity as show on line 13 in the algorithm presented in section 3.1, and also do the same for angular velocity once you solve for the change in orientation?
raigan2 wrote: For constraints related to orientation we define directions as vectors and use atan2, i.e C = atan2(d1) - atan2(d0) would constrain d0 and d1 to be parallel
This is an interesting idea. I also simply measured the angle between the point of contact before and after position correction, which is less work computationally than solving equation (9) with the gradient of C = atan2(d1) - atan2(d0).... Is this not a good method?
raigan2 wrote: Collision is also a pain
I've been using MPR (http://code.google.com/p/mpr2d/), to try something new, which gives me the deepest point of penetration (p), the surface point (q) which is closest to
p and the surface normal (n) at this position. Then I use the constraint, C = (p - q) dot n, as shown in section 3.4, to solve for dP and dQ. This actually gives me very good results when bodies have no angular momentum; I can stack many objects (boxes) with very good stability - visually I detect absolutely no bounce....

However, once I induce angular momentum and solve for change in orientation, and update angular velocity, resting contacts bounce and vibrate all over the place - there doesn't seem to be any loss in kinetic energy. Is this where friction should come into play?
raigan2
Posts: 197
Joined: Sat Aug 19, 2006 11:52 pm

Re: Position Based Dynamics (PBD)

Post by raigan2 »

zzzzrrr wrote:Since the solver is updating position, why also constrain linear velocity? Why not update linear velocity as show on line 13 in the algorithm presented in section 3.1, and also do the same for angular velocity once you solve for the change in orientation?
I simply meant that if you wish to, you can solve at both the position and velocity level in two separate phases, as in Box2D -- this is one possible method to add velocity-level constraints to PBD. This will probably end up being identical to Box2D with NGS position solving.

zzzzrrr wrote: This is an interesting idea. I also simply measured the angle between the point of contact before and after position correction, which is less work computationally than solving equation (9) with the gradient of C = atan2(d1) - atan2(d0).... Is this not a good method?
I don't understand what you mean by "measured the angle between the point of contact"; is this for friction? Do you mean than, rather than decomposing collision constraints into two axial/linear constraints -- one along the surface normal and the other along the tangent -- you're instead using two polar constraints, i.e distance and angle?

zzzzrrr wrote: I've been using MPR (http://code.google.com/p/mpr2d/), to try something new, which gives me the deepest point of penetration (p), the surface point (q) which is closest to
p and the surface normal (n) at this position. Then I use the constraint, C = (p - q) dot n, as shown in section 3.4, to solve for dP and dQ. This actually gives me very good results when bodies have no angular momentum; I can stack many objects (boxes) with very good stability - visually I detect absolutely no bounce....
The main problem with collisions from my experience is that, since position and velocity are coupled, this can add energy -- when you separate penetrating objects you introduce some outward velocity. Since the solver is always moving objects, this means that you can't simply generate collision constraints at the start of the frame (as in Box2D for instance) because later on during solving some new collisions may have been created.

For instance, you'd like your simulation to be something like this:

Code: Select all

Step_Simulation()
{ 
   step world forward in time

   find all colliding pairs A,B

   for(n iterations)
   {
      Solve_Joints()

      for(each pair): Solve_Collision(A,B)
   }
}
However this can cause problems because the solving process can create collisions, which you won't find until next step -- at which points objects are penetrating and pushing them out will add energy to the system.

Something like this will work better:

Code: Select all

Step_Simulation()
{ 
   step world forward in time   

   for(n iterations)
   {
      Solve_Joints()

      for(each body A)
         for(each body B)
            if(A collides with B): Solve_Collision(A,B)
   }
}
will behave better, however:
a) you're now doing a full pass of collision detection at each solver iteration
b) it's still not a perfect solution, because collisions solved later may violate collisions solved earlier (i.e there may still be some residual penetration at the start of the next frame).
c) you can't use any "all-pairs" collision detection method (which are typically faster than pair-by-pair) since you're moving each pair
c) it's not actually possible to do this (as far as I can tell) if you want to avoid O(n^2) collision tests; since A can be moved after the first B, you basically need to re-query your structure every time you move A.. at least, I haven't found a nice way to use broadphase collision systems. Such a way might exist though, and I'm just missing something.


Currently we're doing this:

Code: Select all

Step_Simulation()
{ 
   step world forward in time

   find all *potentially* colliding pairs A,B

   for(n iterations)
   {
      Solve_Joints()

      for(each pair): 
         if(A collides with B): Solve_Collision(A,B)
   }
}
This is basically caching the results of broadphase so that the spatial query system (grid, SAP, etc) doesn't have to be recalculated every time an object is moved. However, the "potentially" colliding pairs (in our current implementation, all objects within some distance d of body A are considered potentially colliding with A) are much more numerous than the actually colliding pairs, and you still need to do narrow-phase collision tests before solving. Also you're still not guaranteed to catch all collisions since the solver may still move A more than d during solving, creating un-predicted collisions.

If you're not seeing problems of this sort, I'd be _very_ interesting in learning specifically how you're processing collisions, since this has so far been the #1 pain of position-based methods.

zzzzrrr wrote: However, once I induce angular momentum and solve for change in orientation, and update angular velocity, resting contacts bounce and vibrate all over the place - there doesn't seem to be any loss in kinetic energy. Is this where friction should come into play?
I'm not sure -- in our case stacks were stable without friction (they would eventually slide/slip apart though), but we're not using proper rigid bodies. I suspect that friction would only help hide whatever the fundamental problem is in this case. Are you modelling bounce/restitution?

We simply decided to forget about bounce and only model totally inelastic collisions to start with, since bounce can always be modelled via gameplay logic (i.e when collision is detected, create a velocity constraint for next frame); in our game the vast majority of objects shouldn't bounce anyway.
zzzzrrr
Posts: 10
Joined: Mon Apr 16, 2007 7:11 pm
Location: Washington, DC

Re: Position Based Dynamics (PBD)

Post by zzzzrrr »

raigan2 wrote:... to constrain the position of two points p0,p1 along an axis parallel to n you would have C = Dot(n,p1 - p0), to constrain the velocity of the points you just use C = Dot(n,v1 - v0). In this case the velocity-level Jacobian is the same as for the position constraints, except that you would replace p with v in the formula.
Here's the velocity constraint from Erin Catto's 2009 GDC slides:
constraint.png
constraint.png (4.69 KiB) Viewed 9255 times
Here's equation (9) from the PBD paper:
PBD.png
PBD.png (4.96 KiB) Viewed 9255 times
Would it be feasible to solve for position, or velocity correction, by using the above constraint equation? From what you are suggesting it would be ok to substitute position for velocity and orientation for angular velocity...?

EDIT:

Using the excellent (and free!) SAGE math tool, I get the following for angular velocity:

C.diff(w2) = (p.y - x2.y)*n.y + (p.x - x2.x)*n.x
C.diff(w1) = -(p.x - x1.y)*n.y - (p.x - x1.x)*n.x

I'm wondering if this would be the same gradient for orientation... e.g. is it safe to assume this:

C.diff(w2) == C.diff(q2)
C.diff(w1) == C.diff(q1)

With q1 and q2 representing orientation.
Dirk Gregorius
Posts: 861
Joined: Sun Jul 03, 2005 4:06 pm
Location: Kirkland, WA

Re: Position Based Dynamics (PBD)

Post by Dirk Gregorius »

You actually end up with a very similar formula for the velocity and the position constraint (using primes instead of dots for simpler reading and typing ):

J*W*JT * lambda = -C
J*W*JT * lambda = -C'

In your example this would be:

C = x2 + r2 - x1 - r1 where r = p - x
C' = J*v = v2 + cross( omega2, r2 ) - v1 - cross( omega1, r1 )

The gotcha is to realize that when solving for the position error and applying the correction (pseudo) impulse doesn't necessarily imply that the constraint is fully satisfied by this since this is a *non-linear* problem. In other words there still can be an error! Remember that you used Taylor expansion:

C(x + dx) ~= C(x) + J * dx

You could take this actually a step further and solve a quadratic formulation (which is unpractical in my opinion, but maybe helps understanding the problem):

C(x + dx) ~= C(x) + J * dx + transpose( dx ) * H * dx

Here H is the Hessian.


Cheers,
-Dirk
Post Reply