Constraints and linear programming

If you have only equality constraints and if these constraints are also all linear independent then solving globally for all constraints can be done in linear time as e.g. shown by D. Baraff.

If you also have inequality constraints like e.g. contacts, limits, or motors and these constraints are also all linear independent then you can extent Baraffs method (though it is not linear time anymore) or use Danzig (e.g. there is a *patented* Danzig solver in the ODE) or Lemke. I think Lemke can also be extended to deal with linear dependent constraints.

If you plan to use you solver in games there is hardly any other option than an iterative solver since you simply cannot guarantee that the current system is solvable. Also global solver are *much* slower than their iterative counterparts and they also need more temporary memory. Finally note that solving an LCP or MLCP is not a linear but a quadratic programming problem.


HTH,
-Dirk

Quote:

Original post by DonDickieD If you have only equality constraints and if these constraints are also all linear independent then solving globally for all constraints can be done in linear time as e.g. shown by D. Baraff.

That sounds familiar. I've probably seen the paper before. Do you have a link for it?

Just from a mathematical standpoint, the matrix for the constraints is extremely sparse, so I was thinking of constructing the matrix using a sparse implementation (basically something that only stores values of the matrix that are non zero), and then solving the system using LU decomposition, but modified so the pivot value is chosen somehow to minimize fill in. Such a system should (as in I think it is the case) be roughly O(n) where n is the number of non zero elements in the matrix.

Since n is related to the number of constraints by (probably) a constant upper bound, an algorithm to solve the system of constraints in linear time would seem possible. I'm curious to compare which methods are more/less numerically stable, memory intensive, etc. so I'm trying to get them all in a common form.

Quote:

If you also have inequality constraints like e.g. contacts, limits, or motors and these constraints are also all linear independent then you can extent Baraffs method (though it is not linear time anymore) or use Danzig (e.g. there is a *patented* Danzig solver in the ODE) or Lemke. I think Lemke can also be extended to deal with linear dependent constraints.

I'm less interested with solvers atm, and more interested in just understanding the underlying math. Like the iterative solvers really are solving an underlying system of equations, so I'd like to play around with things like global solvers, or solving the system using something in between (iterative kernel, I'm thinking). But to do that I need to put the problem into a "familiar" form.

I guess a large part of the problem is that the literature on optimization is pretty dense with confusing notation, and I'm not even sure what to focus on to start learning. So someone telling me that the global system of constraint inequalities can be solved with X optimization technique would at least put me in the right direction to look at literature.

Ideally I'd like to convert the problem into something purely mathematical, and then attempt solving it using techniques presented in the Numerical Recipes book. That is, I want to try and approach this less like a game programmer and more like a mathematician.

Quote:

If you plan to use you solver in games there is hardly any other option than an iterative solver since you simply cannot guarantee that the current system is solvable. Also global solver are *much* slower than their iterative counterparts and they also need more temporary memory.

I'm just goofing around with it because I've never seen a physics engine with a global solver and I think, if nothing else, it's a good learning exercise since it let's me bridge the gap between the purely theoretical nature of the problem and the actual implementation that's commonly done. I'm far less interested in practical implementation at this point.

Just to clarify what I'm after, I have a system of equations and a bunch of box constraints.

Aj = b, where A is a matrix, j is a vector of impulses (in the presentations by eric catto on constraints, these would be the unknown lambda terms), and b is a vector of initial conditions for the system.

I know how to get this far, but now I'm interested in adding two types of constraints:

1. Box constraints. I want to be able to constrain the impulse vector between a "min" vector and a "max" vector like so: min_n <= j_n <= max_n. This needs to be a hard constraint. It's more important to satisfy these constraints than satisfy the Aj = b system.

2. Dependent box constraints. For things like friction, the max/min parts of the box constraint are dependent on another constraint. So basically I have -u * jn <= jf <= u * jn.

And I'm looking for what optimization technique (linear programming, quadratic programming, something more advanced?) I should look into.

I suggest reading all his papers here: http://www.pixar.com/companyinfo/research/deb/

In my above post I was refering to this paper: "Linear-time dynamics using Lagrange multipliers"


If you want to learn the math I recommend these books:

http://www.amazon.com/Solving-Ordinary-Differential-Equations-1/dp/B001JPQ4ZK/ref=sr_1_6?ie=UTF8&s=books&qid=1247302060&sr=1-6

http://www.amazon.com/Solving-Ordinary-Differential-Equations-Differential-Algebraic/dp/0387537759/ref=sr_1_3?ie=UTF8&s=books&qid=1247302060&sr=1-3

http://www.amazon.com/Geometric-Numerical-Integration-Structure-Preserving-Computational/dp/3540306633/ref=pd_sim_b_2

I assume you want to learn this in the context of rigid body dynamics. So you have to solve (ignoring inequality for simpliicty right now);

dv/dt = ( f_ext + JT * lambda ) / M
dx/dt = v

C(x) = 0

The problem is that the constraint forces (lambda) are not known a priori. There exist several possibilites to solve this problem. You can solve on the acceleration level like Brarff does, but you can suffer from infinite friction forces. Baraff claims that this is not a real problem and I am also not convinced that you run into problems. This approach transforms the DAE problem into an ODE and you can use whatever integrator you like since you solve for the constraint forces explicetly. The other approach is to solve on the velocity level. The was e.g. introduced by Anitescu and is used in all engines I am aware of (e.g. ODE, Bullet, Box2D and also the commercial ones). In the sequential impulse implementation you basically filter out the velocites that were otherwise canceled out by the constraint forces. Finally you can solve on the position level like e.g. introduced by T. Jacobsen and M. Mueller in his paper "Position Based Dynamics". You can find it e.g. here:

http://www.matthiasmueller.info/

Of course things get more complicated when you add inequatlity constraints. Gino v. d. Bergen and myself are editors on a book about game physic pearls and this topic will be covered there in detail by authors like Erin Catto, Thomas Jacobsen and myself. Hopefully can bring some light into this topic.


HTH,
-Dirk

Quote:

Original post by DonDickieD I suggest reading all his papers here: http://www.pixar.com/companyinfo/research/deb/ In my above post I was refering to this paper: "Linear-time dynamics using Lagrange multipliers"

Thanks, I'll give it a look.

Quote:

If you want to learn the math I recommend these books: http://www.amazon.com/Solving-Ordinary-Differential-Equations-1/dp/B001JPQ4ZK/ref=sr_1_6?ie=UTF8&s=books&qid=1247302060&sr=1-6 http://www.amazon.com/Solving-Ordinary-Differential-Equations-Differential-Algebraic/dp/0387537759/ref=sr_1_3?ie=UTF8&s=books&qid=1247302060&sr=1-3 http://www.amazon.com/Geometric-Numerical-Integration-Structure-Preserving-Computational/dp/3540306633/ref=pd_sim_b_2 I assume you want to learn this in the context of rigid body dynamics. So you have to solve (ignoring inequality for simpliicty right now); dv/dt = ( f_ext + JT * lambda ) / M dx/dt = v C(x) = 0 The problem is that the constraint forces (lambda) are not known a priori. There exist several possibilites to solve this problem. You can solve on the acceleration level like Brarff does, but you can suffer from infinite friction forces. Baraff claims that this is not a real problem and I am also not convinced that you run into problems. This approach transforms the DAE problem into an ODE and you can use whatever integrator you like since you solve for the constraint forces explicetly. The other approach is to solve on the velocity level. The was e.g. introduced by Anitescu and is used in all engines I am aware of (e.g. ODE, Bullet, Box2D and also the commercial ones). In the sequential impulse implementation you basically filter out the velocites that were otherwise canceled out by the constraint forces. Finally you can solve on the position level like e.g. introduced by T. Jacobsen and M. Mueller in his paper "Position Based Dynamics". You can find it e.g. here: http://www.matthiasmueller.info/ Of course things get more complicated when you add inequatlity constraints. Gino v. d. Bergen and myself are editors on a book about game physic pearls and this topic will be covered there in detail by authors like Erin Catto, Thomas Jacobsen and myself. Hopefully can bring some light into this topic.

Just to be clear, I took linear algebra, ODE, and numerical analysis classes in college. I know how to formulate and solve equations of the form y''(t) := a(t) + b * y'(t) + k * y(t) for instance, both numerically and analytically. And I know how to form and solve systems of linear equations. So I do have a pretty solid grounding. But I never learned anything about optimization beyond optimizing single dimensional functions in calculus. So there's this giant gaping hole in my knowledge.

As for constraints, I've worked through Erin Catto's presentations for GDC. So I understand the idea of formulating the constraint function, finding it's derivative, etc. And I understand setting up the Jacobian (both the local Jacobian and the global Jacobian) and solving for the unknown scalar impulse terms (lambda).

I'm doing a simple 2D physics engine for work right now (as a training excercise). Right now the only constraints are frictionless contact constraints with baumgarte correction. By solving the contact manifolds of each body pair at once (instead of sequentially), I've gotten really impressive results. Since it's 2D, the contact manifold forms a 2x2 system of equations (for stacking blocks anyway), and I solve for the two impulse terms using Cramer's rule. I've managed to get stable stacking for upwards of 100 bodies with no visible overlap, and when it does become unstable, it's from numerical errors accumulating laterally (not vertically), and since there's no friction to slow them down they pop out like wet soap. Even still, the stack is stable for several seconds.

I could naively add friction terms as a second system of equations, and treat them differently, but if you do the analysis friction terms and contact impulse terms aren't independent. Friction can influence contact impulses beyond even the constraints on friction. I figure I'll get more stable results by solving the friction at the same time as the contact impulses. If I solve for friction as if there's no constraint on its value, it's still very straightforward and I have a 4x4 system to solve. Still quite manageable.

The trick now becomes adding in the constraints. Making sure that contacts only push, not pull (non negativity). And making sure that friction is bounded by its associated contact impulse. I'm not sure, but I'm assuming you can't just apply the constraints as an afterthought without getting potentially incorrect results.

I'm sure there's an optimization technique which applies, but I really don't know it. I'm trying to make my way through optimization literature in the Numerical Recipes book and another book from work, but the literature is dense. And linear programming, at least, seems to assume that the unknown vector your solving for should have nonnegativity, which isn't something I want to enforce for all constraints. Especially not friction.

And of course once I figure out just the case of just two rigid bodies colliding, I should be able to expand it to solve the entire world at once as a single optimization problem. A sparse-smart solver should be able to solve the system in number of constraints linear time, I think.

In the presence of unilateral constraints you have to solve an LCP instead of a LS. The LCP can be transformed into a quadratic optimization problem. IIRC they are equivalent. Friction is an additional problem since the friction constraints depent on the non-penetration constraint and make the effective mass matrix kind of ill-behaved. Look at Murty's webbook here:

http://ioe.engin.umich.edu/people/fac/books/murty/linear_complementarity_webbook/

I would concentrate more on DAE and LCP instead on optimization, but this is my very personal opinion on this. If you are familiar with Box2D look at e.g. the contact block solver or the revolute joint with limits. These are both interesting implementations of the "Direct Enumeration Method" as explained in the book above at the very beginning.

If you want to understand optimization you first have to re-formulate the problem in my opinion. I think this is possible, but I have never done this.


Cheers,
-Dirk

Quote:

Original post by DonDickieD In the presence of unilateral constraints you have to solve an LCP instead of a LS. The LCP can be transformed into a quadratic optimization problem. IIRC they are equivalent. Friction is an additional problem since the friction constraints depent on the non-penetration constraint and make the effective mass matrix kind of ill-behaved. Look at Murty's webbook here: http://ioe.engin.umich.edu/people/fac/books/murty/linear_complementarity_webbook/

Thanks, I'll give it a read, too.

Quote:

I would concentrate more on DAE and LCP instead on optimization, but this is my very personal opinion on this. If you are familiar with Box2D look at e.g. the contact block solver or the revolute joint with limits. These are both interesting implementations of the "Direct Enumeration Method" as explained in the book above at the very beginning. If you want to understand optimization you first have to re-formulate the problem in my opinion. I think this is possible, but I have never done this.

Okay, let's assume for the moment I want to formulate the constraints using LCP and solve it that way. Suppose I have a pin constraint in 2D, so it's basically like a contact constraint but it can pull as well as push. And suppose I want to add limits to that pin constraint so that it's impulse (lambda) term is bounded by [-100, 50].

Suppose it's the only constraint in the world. How do you convert it in to a form LCP can solve? Because as near as I can tell LCP has that non negativity constraint built in just like linear programming does.

If you could walk through the steps of solving just this single constraint using LCP I can probably take it from there (that is, assume that the M matrix is still a matrix even though it has only a single element, etc.)

Solving a 1D constraint in rigid body dynamics is straigt forward. Actually this is what you do all the time in the iterative solvers. You have

w = a * x - b >= 0 and x >= 0 and w * x = 0

In rigid body dynamics the effective mass matrix (J*W*JT)^-1 is PD (and SPD in the presence of friction). This means that a > 0 (Practical hint: You should assert this in your code!)

So you have x = b / a (with a > 0 as explained above)

Case 1 (b >= 0): In this case x >= 0 and therefore w = 0
Case 2 (b < 0): In this case x < 0 and therefore we clamp x = 0 -> w > 0


You have examples directly in the beginning of the web book I mentioned. Also Box2D has the block-solver in the contact code and revolute joint with limits. Note that Erin shifts variables to account for the clamping against the accumulated impulse which is a small additional gotcha.


HTH,
-Dirk

Don't remember how to get from an LCP to a boxed LCP right now. I think you have to shift variables. Maybe google for box LCP.

Also make sure that you understand what w = a * x - b >= 0 and x >= 0 and w * x = 0 means physically. Can you explain what w and x are? This helps alot for the understanding.


HTH,
-Dirk

Quote:

Original post by Numsgil I solve for the two impulse terms using Cramer's rule.

FYI, Cramer's rule is a particularly slow way to solve linear systems of equations. Gaussian elimination with pivoting is the standard non-iterative method.

Quote:

And linear programming, at least, seems to assume that the unknown vector your solving for should have nonnegativity, which isn't something I want to enforce for all constraints. Especially not friction.

I don't know much about simulating collisions, but I can comment on formulating LPs (linear programming problems). There are some very simple ways to deal with this.

If a decision variable x can be either positive or negative it is called unrestricted, and you can express it as x = x₊ - x_- where x₊ and x_- are both nonnegative. Alternatively, if there is a lower bound (call it L) on x then you can of course replace x by L + x_d where x_d is the nonnegative amount by which x exceeds its lower bound.

I admit that these both feel a little hackish, but they're common practice.

Hmmm; I hadn't realized this before: Actually, there's an even better way to handle unrestricted variables when there's no lower bound, which doesn't make the dimension of your problem increase so much. You only need to add one decision variable.

Say you have a homogeneous LP,

min a x + b y + c z

st

d x + e y + f z > 0
g x + h y + i z > 0 .

where x,y,z are your decision variables and a,b,c,d,e,f,g,h,i are just scalar constants.

It's called homogeneous because the constraints are all of the form "[linear combination of decision variables] > 0" instead of "[linear combination of decision variables] > nonzero." This is general since you can put any problem into homogeneous form by introducing slack variables.

Now, express your decision variables as

x = X - L
y = Y - L
z = Z - L

where X,Y,Z,L are all nonnegative decision variables (so we've introduced one new variable, L). Then the minimization becomes

min a (X - L) + b (Y - L) + c (Z - L) =
min a X + b Y + c Z - (a + b + c)L .

As for our constraints...

The first constraint,
d (X - L) + e (Y - L) + f (Z - L) > 0
reduces to
d X + e Y + f Z - (d + e + f) L > 0

and the second constraint,
g (X - L) + h (Y - L) + i (Z - L) > 0
reduces to
g X + h Y + i Z - (g + h + i) L > 0

so essentially you just need to add up the elements of each row of your constraint matrix, flip the sign, and append these numbers as an extra column.

In this way you can handle as many unrestricted variables as you like while only introducing one new variable.

Quote:

Original post by Emergent Quote: Original post by Numsgil I solve for the two impulse terms using Cramer's rule. FYI, Cramer's rule is a particularly slow way to solve linear systems of equations. Gaussian elimination with pivoting is the standard non-iterative method.

You mean that Cramer's rule is a particularly non scalable way of solving linear systems, and shouldn't be used for more than about a 4x4 system, but it's a particularly fine choice for a 2x2 system, and that LU decomposition with partial pivoting is the standard non-iterative method (since it requires a whopping 3 times fewer operations for any given system size), right? ;)

Quote:

If a decision variable x can be either positive or negative it is called unrestricted, and you can express it as x = x+ - x- where x+ and x- are both nonnegative. Alternatively, if there is a lower bound (call it L) on x then you can of course replace x by L + xd where xd is the nonnegative amount by which x exceeds its lower bound. I admit that these both feel a little hackish, but they're common practice.

Huh, that's way simpler than I would have guessed. Thanks, that's at least one major hurdle to my understanding out of the way.

Quote:

Original post by Numsgil You mean that Cramer's rule is a particularly non scalable way of solving linear systems, and shouldn't be used for more than about a 4x4 system, but it's a particularly fine choice for a 2x2 system, and that LU decomposition with partial pivoting is the standard non-iterative method (since it requires a whopping 3 times fewer operations for any given system size), right? ;)

Ok, cool; you're up on this. Nice! :-)

Since you'd mentioned some basic calculus experience I've got to admit that I was picturing a straightforward implementation of the algorithms that you get taught in a course like that: E.g., "Compute the inverse using Cramer's rule, then multiply by the inverse -- and since Cramer's rule requires the computation of determinants, compute them recursively using Laplace's formula." That's pretty much the worst way to solve equations in the world, but I think it's what I'd expect somebody to do if they didn't know better.

And actually, I should admit that I'm not all that clear on LU decomposition "vs." Gaussian elimination myself; I don't "get" the distinction that is being made. My understanding is that Gaussian elimination just does the operations in place, whereas LU decomposition constructs upper and lower triangular matrices that "do book-keeping" for what occurs during Gaussian elimination. I understand the advantage of computing LU when you need to solve systems with the same coefficients over and over (actually it's better to have L and U than the inverse; I get this), but if you only need to solve a system once it'd seem to me that plain old Gaussian elimination (with (partial) pivoting) would be the better choice since it does everything in place. Have I missed the boat on something?

Quote:

Huh, that's way simpler than I would have guessed. Thanks, that's at least one major hurdle to my understanding out of the way.

Glad I was able to tell you something new!

Quote:

Original post by Emergent And actually, I should admit that I'm not all that clear on LU decomposition "vs." Gaussian elimination myself; I don't "get" the distinction that is being made. My understanding is that Gaussian elimination just does the operations in place, whereas LU decomposition constructs upper and lower triangular matrices that "do book-keeping" for what occurs during Gaussian elimination. I understand the advantage of computing LU when you need to solve systems with the same coefficients over and over (actually it's better to have L and U than the inverse; I get this), but if you only need to solve a system once it'd seem to me that plain old Gaussian elimination (with (partial) pivoting) would be the better choice since it does everything in place. Have I missed the boat on something?

The first trick is that you can do the LU decomposition in place same as Guassian Elimination, so there isn't any extra memory overhead for LU decomposition. The second trick is that LU decomposition means you don't have to compute the inverse of the matrix to solve the system. Gaussian Elimination, from my understanding, solves for the inverse of the matrix and then uses it to solve for the unknowns. But using LU decomposition you only have to decompose the matrix, which saves computational time.

If you actually need the real inverse of the matrix, both methods are essentially the same. They're both O(n^3) operations (no free lunch), but it just so happens that if you just want to solve a system of equations, LU decomposition is a constant order of magnitude faster (3x faster) than finding the inverse.

There are other methods that can give better (faster) results for certain problem types, but the workhorse for solving linear systems is LU decomposition. As for the pivoting, partial pivoting is almost as good as full pivoting, but is much easier to actually build and use.

Chapter 2 of Numerical Recipes in C goes over all of this (plus way more I haven't even looked at yet) and provides some sample implementations.

....

I'm still working my way through the optimization literature. I mostly understand Linear Programming and the different methods for solving it. I haven't really looked at LCP, but when I do I'll post any questions I have about it.

[Edited by - Numsgil on July 13, 2009 11:19:16 AM]

Quote:

Original post by DonDickieD The other approach is to solve on the velocity level. The was e.g. introduced by Anitescu and is used in all engines I am aware of (e.g. ODE, Bullet, Box2D and also the commercial ones).

Do you have any references to papers by Anitescu on this topic, say, ones that might summarize the ideas?

Quote:

Original post by DonDickieD Solving a 1D constraint in rigid body dynamics is straigt forward. Actually this is what you do all the time in the iterative solvers. You have w = a * x - b >= 0 and x >= 0 and w * x = 0 In rigid body dynamics the effective mass matrix (J*W*JT)^-1 is PD (and SPD in the presence of friction). This means that a > 0 (Practical hint: You should assert this in your code!) So you have x = b / a (with a > 0 as explained above) Case 1 (b >= 0): In this case x >= 0 and therefore w = 0 Case 2 (b < 0): In this case x < 0 and therefore we clamp x = 0 -> w > 0 You have examples directly in the beginning of the web book I mentioned. Also Box2D has the block-solver in the contact code and revolute joint with limits. Note that Erin shifts variables to account for the clamping against the accumulated impulse which is a small additional gotcha. HTH, -Dirk

I was hoping to make my way through your explanation, but you're sort of cheating by using different notation than the book you directed me to, and assuming that your matrix is a single scalar, which obviously won't work as the system size increases.

So let me try to jive what you wrote and my understanding of LCP a bit:

Suppose I can represent the unconstrained system as a system of equalities like so:

| O11 O12 | * |j0| = |b0|| O21 O22 |   |j1|   |b1|

That is, it's a 2x2 system. I can write this more succinctly as: O * j = b. (O,b) are known, j is the impulse vector I want to find.

Back to the LCP, the canonical form is: given a matrix M and a vector q, find w, z such that:

w - Mz = q
w,z non-negative and complimentary

First question: how do (O,b) and (M,q) relate? From what I can tell, they would be the same, except that the LCP would now impose a non negativity constraint on the solution.

Second question: Likewise, is z the same thing as my j vector? From what I can tell w is essentially a slack vector and you would just ignore it in extracting the impulses you want from z.

Let's consider the slightly more interesting case of two bodies in a chain in collision with a static ceiling:

          Ceiling______________________________            |   j0       ------------       |          |                     |  Body A  |   vel A       |          |     ||       ------------     \/            ||  j1       ------------       |          |       |  Body B  |   vel B       |          |     ||       ------------     \/

Just by looking at it empirically, we can tell that j0 and j1 should be 0. Let's play with some nice numbers and let the system be:

| 1 1 | * |j0| = |3|| 1 2 |   |j1|   |4|

Solving with Cramer's rule we get j0 = 2, j1 = 1 (where + is assumed to be towards the ceiling). This isn't the correct answer, of course, because we haven't added the condition that j0 <= 0 yet.

So how do I convert this system in to something LCP can solve?

@Numsgil:
My example explains how to solve 1D constraint to get you the idea. I posted several links which give you other approaches how to solve for more than 1 constraint, e.g. Direct Enumeration, Danzig, Lemke, and Baraff.

@Dave
Look at Erin Catto's first presentation at the GDC and e.g. Anitescu PhD (I think it is not freely available - contact me if you cannot find it). Erin's presentation can be found here:

http://www.gphysics.com/files/IterativeDynamics.pdf

Quote:

Original post by DonDickieD @Numsgil: My example explains how to solve 1D constraint to get you the idea. I posted several links which give you other approaches how to solve for more than 1 constraint, e.g. Direct Enumeration, Danzig, Lemke, and Baraff.

But that's how to solve the problem once it's in the LCP format. I'm still trying to understand how to format it that way in the first place. I'd really appreciate you going over my post and at least answering the first two questions about how the w - Mz = q format corresponds to solving the problem without constraints.

I've tried googling, but since most real time engines use basically a sequential solver any relationship to LCP is more a passing reference than a formal exploration. And the few formal explorations I found (siggraph papers, etc.) were extremely dense with foreign notation. I appreciate the links to webbooks, etc. you've posted, and I'm still going through them all, but an actual concrete example would go a long way towards just giving me the big picture, which I'm still missing.

I'm using a slightly different form for the example I gave, so feel free to reformulate it however you like. The chain in collision contact with the ceiling is still a pretty illustrative example and it would answer a lot of my questions to see it solved using LCP.

Quote:

Original post by DonDickieD Look at Erin Catto's first presentation at the GDC and e.g. Anitescu PhD (I think it is not freely available - contact me if you cannot find it). Erin's presentation can be found here: http://www.gphysics.com/files/IterativeDynamics.pdf

Thanks, I already have Catto's PDF. I have looked a long time at Anitescu's web site, but cannot see anything labeled as his PhD thesis (the university page for graduated students has no links). I do not mind paying to download the thesis, but I have to find a link and a way to do this first :)