Modeling and Solving Constraints. Basic Idea
1. Modeling and Solving ConstraintsErin Catto
2. Basic IdeaConstraints are used to simulate joints, contact,
We need to solve the constraints to stack boxes
and to keep ragdoll limbs attached.
Constraint solvers do this by calculating impulse
or forces, and applying them to the constrained
3. OverviewConstraint Formulas
Jacobians, Lagrange Multipliers
Joints, Motors, Contact
Building a Constraint Solver
4. Constraint TypesContact and Friction
5. Constraint TypesRagdolls
6. Constraint TypesParticles and Cloth
7. Show Me the Demo!
8. Bead on a 2D Rigid WireImplicit Curve Equation:
This is the position constraint.
C ( x, y ) 0
9. How does it move?The normal vector is perpendicular to the velocity.
dot(n, v) 0
10. Enter The CalculusPosition Constraint:
C ( x) 0
If C is zero, then its time derivative is zero.
11. Velocity ConstraintC 0
Velocity constraints define the allowed motion.
Next we’ll show that velocity constraints depend
linearly on velocity.
12. The JacobianDue to the chain rule the velocity constraint has
a special structure:
J is a row vector called the Jacobian.
J depends on position.
The velocity constraint is linear.
13. The JacobianThe Jacobian is perpendicular to the velocity.
C Jv 0
14. Constraint ForceAssume the wire is frictionless.
What is the force between the wire and the bead?
15. Lagrange MultiplierIntuitively the constraint force Fc is parallel to
the normal vector.
16. Lagrange MultiplierThe Lagrange Multiplier (lambda) is the constraint force
We use a constraint solver to compute lambda.
17. Jacobian as a CoordinateTransformSimilar to a rotation matrix.
Except it is missing a couple rows.
So it projects some dimensions to zero.
The transpose is missing some columns, so some
dimensions get added.
18. Velocity Transformv
19. Force TransformJ
20. Refresher: Work and PowerWork = Force times Distance
Work has units of Energy (Joules)
Power = Force times Velocity (Watts)
P dot F, V
21. Principle of Virtual WorkPrinciple: constraint forces do no work.
We can ensure this by using:
Proof (compute the power):
Pc F v J v Jv 0
The power is zero, so the constraint does no work.
22. Constraint QuantitiesPosition Constraint
23. Why all the Painful Abstraction?We want to put all constraints into a common form for the
This allows us to efficiently try different solution
24. Addendum: Modeling Time DependenceSome constraints, like motors, have prescribed motion.
This is represented by time dependence.
C x, t 0
C Jv b(t ) 0
25. Example: Distance Constraintx
is the tension
C x L
26. Gory DetailsdC d
x2 y 2 L
2 x y
2 xvx yv y
2 x y
x v x xT
x y y x
27. Computing the JacobianAt first, it is not easy to compute the Jacobian.
It gets easier with practice.
If you can define a position constraint, you can find its
Here’s how …
28. A Recipe for JUse geometry to write C.
Differentiate C with respect to time.
Identify J and b by inspection.
C Jv b
29. Constraint PotpourriJoints
30. Joint: Distance Constraintx
Fc J T
31. MotorsA motor is a constraint with limited force (torque).
C sin t
Note: this constraint does work.
32. Velocity Only MotorsExample
Usage: A wheel that spins at a constant rate.
We don’t care about the angle.
33. Inequality ConstraintsSo far we’ve looked at equality constraints (because
they are simpler).
Inequality constraints are needed for contact and joint
We put all inequality position constraints into this form:
C (x, t ) 0
34. Inequality ConstraintsThe corresponding velocity constraint:
35. Inequality ConstraintsForce Limits:
Inequality constraints don’t suck.
36. Contact ConstraintNon-penetration.
Friction: sliding, sticking, and rolling
37. Non-Penetration Constraintbody 2
38. Non-Penetration ConstraintC ( v p 2 v p1 ) n
v 2 ω 2 p x 2 v1 ω1 p x1 n
p x n
A B C
C A B
B C A
39. RestitutionRelative normal velocity
( v p 2 v p1 ) n
Adding bounce as a velocity bias
C v ev 0
40. Friction ConstraintFriction is like a velocity-only motor.
The target velocity is zero.
C vp t
v ω p x t
41. Friction ConstraintThe friction force is limited by the normal force.
n t n
3D is a bit more complicated. See the references.
42. Constraints SolversWe have a bunch of constraints.
We have unknown constraint forces.
We need to solve for these constraint forces.
There are many ways different ways to compute
43. Constraint Solver TypesGlobal Solvers (slow)
Iterative Solvers (fast)
44. Solving a Chain1
solve for 1, 2, and 3
solve for 1
solve for 2
solve for 3
45. Sequential Impulses (SI)An iterative solver.
SI applies impulses at each constraint to correct the
SI is fast and stable.
Converges to a global solution.
46. Why Impulses?Easier to deal with friction and collision.
Lets us work with velocity rather than acceleration.
Given the time step, impulse and force are
47. Sequential ImpulsesStep1:
Integrate applied forces, yielding tentative
Apply impulses sequentially for all constraints,
to correct the velocity errors.
Use the new velocities to update the positions.
48. Step 1: Newton’s LawWe separate applied forces and
Mv Fa Fc
49. Step 1: Mass MatrixParticle
m 0 0
M 0 m 0
0 0 m
May involve multiple particles/bodies.
50. Step 1: Applied ForcesApplied forces are computed according to some law.
Gravity: F = mg
Spring: F = -kx
Air resistance: F = -cv2
51. Step 1 : Integrate Applied ForcesEuler’s Method for all bodies.
v 2 v1 hM Fa
This new velocity tends to violate the velocity
52. Step 2: Constraint ImpulseThe constraint impulse is just the time step
times the constraint force.
53. Step 2: Impulse-MomentumNewton’s Law for impulses:
M v Pc
In other words:
v 2 v 2 M Pc
54. Step 2: Computing LambdaFor each constraint, solve these for :
v 2 v 2 M 1Pc
Jv 2 b 0
Note: this usually involves one or two bodies.
55. Step 2: Impulse SolutionmC Jv 2 b
The scalar mC is the effective mass seen by
the constraint impulse:
56. Step 2: Velocity UpdateNow that we solved for lambda, we can use it
to update the velocity.
v 2 v 2 M Pc
Remember: this usually involves one or two bodies.
57. Step 2: IterationLoop over all constraints until you are done:
- Fixed number of iterations.
- Corrective impulses become small.
- Velocity errors become small.
58. Step 3: Integrate PositionsUse the new velocity to integrate all
body positions (and orientations):
x2 x1 hv 2
This is the symplectic Euler integrator.
59. Extensions to Step 2Handle position drift.
Handle force limits.
Handle inequality constraints.
60. Handling Position DriftVelocity constraints are not obeyed precisely.
Joints will fall apart.
61. Baumgarte StabilizationFeed the position error back into the velocity constraint.
New velocity constraint:
62. Baumgarte StabilizationWhat is the solution to this?
First-order differential equation …
C C0 exp
64. Tuning the Bias FactorIf your simulation has instabilities, set the bias factor to
zero and check the stability.
Increase the bias factor slowly until the simulation
Use half of that value.
65. Handling Force LimitsFirst, convert force limits to impulse limits.
impulse h force
66. Handling Impulse LimitsClamping corrective impulses:
clamp , min , max
Is it really that simple?
67. How to ClampEach iteration computes corrective impulses.
Clamping corrective impulses is wrong!
You should clamp the total impulse applied over the
The following example shows why.
68. Example: 2D Inelastic Collisionv
A Falling Box
69. Iterative Solutioniteration 1
Suppose the corrective impulses are too strong.
What should the second iteration look like?
70. Iterative Solutioniteration 2
To keep the box from bouncing, we need
downward corrective impulses.
In other words, the corrective impulses are
71. Iterative SolutionBut clamping the negative corrective impulses
wipes them out:
clamp( , 0, )
This is one way to introduce jitter into
72. Accumulated ImpulsesFor each constraint, keep track of the total impulse
This is the accumulated impulse.
Clamp the accumulated impulse.
This allows the corrective impulse to be negative yet the
accumulated impulse is still positive.
73. New Clamping Procedure1. Compute the corrective impulse, but don’t
2. Make a copy of the old accumulated impulse.
3. Add the corrective impulse to the accumulated
4. Clamp the accumulated impulse.
5. Compute the change in the accumulated
impulse using the copy from step 2.
6. Apply the impulse delta found in Step 5.
74. Handling Inequality ConstraintsBefore iterations, determine if the inequality constraint is
If it is inactive, then ignore it.
Clamp accumulated impulses:
75. Inequality ConstraintsA problem:
Aiming for zero overlap leads to JITTER!
76. Preventing OvershootAllow a little bit of penetration (slop).
If separation < slop
C Jv slop
Note: the slop will be negative (separation).
77. Warm StartingIterative solvers use an initial guess for the lambdas.
So save the lambdas from the previous time step.
Use the stored lambdas as the initial guess for the new
Benefit: improved stacking.
78. Step 1.5Apply the stored impulses.
Use the stored impulses to initialize the accumulated
79. Step 2.5Store the accumulated impulses.
80. Further Reading & Sample CodeFurther Reading &
81. Box2DAn open source 2D physics engine.
Written in C++.