Constraints and linear programming

This topic is 3098 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

Recommended Posts

I'm slowly making my way through the math for constraints, and trying to build a global solver instead of the iterative one most engines use. What's the best way mathematically to solve the constraints all at once? Suppose I've set up my system of equations like this:
| O11 O12 O12 ... |   | j0 |   | b0 |
| O21 O22 O23 ... | * | j1 | = | b1 |
| O31 O32 O33 ... |   | j2 |   | b2 |
|      ...        |   | ...|   | ...|

The jn are basically the unknown scalar impulse of some particular constraint (It's formulated differently from how it's usually done with Jacobians, etc., but I think it's mathematically equivalent). I can solve that system right now and get the impulses for my constraints. But now I want to add things like friction (secondary impulses) and max/min box constraints. Which sounds suspiciously like some form of linear (or quadratic, etc.) programming. So basically I want to add additional box constraints such as: c0.min <= j0 <= c0.max And also add in secondary impulses (the friction) like so: jf <= j0 * coefficient friction plus some sort of constraint saying that friction can't make an object slide backwards. My current thinking is to first solve the system I present at the top to get all the impulse terms, and then set up a linear programming problem to minimize the difference between the calculated constraints and what they need to be to satisfy the additional box constraints. But that sounds really inefficient since I'm basically solving two slightly different linear systems. Plus I really don't know how to formulate friction mathematically so it can fit in to the system. My thinking is that there's a more complex form of linear programming that will do what I want, and I'm just not familiar with it. Can anyone point me in the right direction for this stuff? I've seen Chris Hecker's presentation on MLCP, and my guess is that there's something here that would help me, but I'm having trouble making heads or tales of it. I also have the Numerical Recipes book, and have been trying to make my way through the linear programming sections. If I could even figure out the optimal way to present the system to be solved, that'd get me a long way.

Share on other sites
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

Share on other sites
Quote:
 Original post by DonDickieDIf 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.

Share on other sites
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.

Share on other sites
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

Share on other sites
Quote:
 Original post by DonDickieDI 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-6http://www.amazon.com/Solving-Ordinary-Differential-Equations-Differential-Algebraic/dp/0387537759/ref=sr_1_3?ie=UTF8&s=books&qid=1247302060&sr=1-3http://www.amazon.com/Geometric-Numerical-Integration-Structure-Preserving-Computational/dp/3540306633/ref=pd_sim_b_2I 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) = 0The 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.

Share on other sites
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

Share on other sites
Quote:
 Original post by DonDickieDIn 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.)

Share on other sites
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

Share on other sites
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

Share on other sites
Quote:
 Original post by NumsgilI 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 + 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.

Share on other sites
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.

Share on other sites
Quote:
Original post by Emergent
Quote:
 Original post by NumsgilI 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.

Share on other sites
Quote:
 Original post by NumsgilYou 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!

Share on other sites
Quote:
 Original post by EmergentAnd 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]

Share on other sites
Quote:
 Original post by DonDickieDThe 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?

Share on other sites
Quote:
 Original post by DonDickieDSolving a 1D constraint in rigid body dynamics is straigt forward. Actually this is what you do all the time in the iterative solvers. You havew = a * x - b >= 0 and x >= 0 and w * x = 0In 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 = 0Case 2 (b < 0): In this case x < 0 and therefore we clamp x = 0 -> w > 0You 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?

Share on other sites
@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

Share on other sites
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.

Share on other sites
Quote:
 Original post by DonDickieDLook 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 :)

Share on other sites
Quote:
Original post by Dave Eberly
Quote:
 Original post by DonDickieDLook 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 :)

This page has Anitescu's 1997 University of Iowa PhD thesis as "Multi-Rigid Body Dynamics with Contact and Friction" under Florian Potra.

So, I'm just guessing here, but probably any of Anitescu's paper with a similar title will contain the same results. That is, after all, how the "academic method" works: you publish the same thing in as many journals as possible, under a similar (but not quite identical) name, to get as much publishing credit as possible.

If your search for that title, you find e.g. this ppt by Anitescu which references 'Anitescu & Potra, 1996' for "Newton's law with impulses and velocities" (as well as Stewart & Trinkle, 1995).

Guessing again, that 1996 reference is either one of these (both an off-by-one error to that reference):

38. M. Anitescu and F.A.Potra, Formulating dynamic multi-rigid-body contact problems with friction as solvable Linear Complementarity Problems Nonlinear Dynamics, 14, pp. 231–247, 1997.
39. M. Anitescu, J. F. Cremer and F. A. Potra, Formulating 3D Contact Dynamics Problems, Mechanics of Structures and Machines, 24(4), November 1996.

Share on other sites
Assume you have a distance constraint between two particles. Also assume that we only can push the particles apart.

C(x1, x2) = |x2 - x1| - L = 0
dC/dt = (x2 - x1) / |x2 - x1| * (v2 - v1)

Define n = (x2 - x1) / |x2 - x1| -> dC/dt = n * (v2 - v1)

We apply impulses P1 and P2:

(1) v1' = v1 + P1 / m1
(2) v2' = v2 + P2 / m2

such that the post-velocities v1' and v2' should satisfy the velocity constraint dC/dt:

(3) n * (v2' - v1') = 0

We also know the impulse direction (results from the d'Alembert Principle):

(4) P1 = -lambda * n
(5) P2 = lambda * n

5 equations for 5 unknowns(v1', v2', P1, P2, lambda). Now we solve for lambda using a more compact vector notation:

v' = v + P / M
J * v' = 0
P = JT * lambda

<=> J * M^-1 * JT * lambda = -J * v

The above equation finds an impulse lambda such that the relative velocites along n are cancelled out. Now let's form an LCP from this:

v_rel' = J * M^-1 * JT * lambda + J * v >= 0 _and_
lambda >= 0 _and_
v_rel' * lambda = 0

If the particles are drifting apart lambda becomes less than zero and we ignore the impulse. If the particles are approaching each other we apply a positive impulse lambda such that v_rel' becomes zero.

v_rel' is the relative velocity after we apply the impulse. J * v is the relative velocity before we apply the impulse. If the relative velocity is greater zero lambda becomes zero. Otherwise we compute a lambda such that the relative velocity after applying the impulse becomes zero.

Note that LCPs and sequential impulses are mathemetical equivalent. It is basically just the physical interpretation what a projected Gauss-Seidel does when solving a MLCP. Namely it applies "impulses sequentially".

HTH,
-Dirk

Share on other sites
Quote:
 Original post by Christer EricsonThis page has Anitescu's 1997 University of Iowa PhD thesis as "Multi-Rigid Body Dynamics with Contact and Friction" under Florian Potra.

Thanks. I had found the reference, just not a link to a copy of it. The PPT file at the link you posted is, well, a PPT file--not much detail. Usually, a thesis has a lot more details than you would find in a research publication or presentation, so I was hoping to read the original.

Quote:
 So, I'm just guessing here, but probably any of Anitescu's paper with a similar title will contain the same results. That is, after all, how the "academic method" works: you publish the same thing in as many journals as possible, under a similar (but not quite identical) name, to get as much publishing credit as possible.

We called this practice "splitsmanship" :) I never practiced it and I don't miss the publish-or-perish environment. I'll track down one of Anitescu's earlier papers and see what I can find.

Share on other sites
Quote:
 Original post by NumsgilGaussian Elimination, from my understanding, solves for the inverse of the matrix and then uses it to solve for the unknowns.

Thanks for your reply, Numsgil. But FYI -- and, apologies for this somewhat off-topic parenthetical discussion in this thread -- this is not correct. Solving a single system of equations with Gaussian elimination does not involve constructing an inverse. Just the opposite: You compute the inverse with GE by solving n equations (one for each element of the natural basis).

Quote:
 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.

AFAIK, using LU decomposition to compute an nxn inverse is O(n^3). Using Gaussian elimination to do so is O(n^4).

Quote:
 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.

Thanks for the reference. The Numerical Recipes people make me sad though: The books used to be available for free online as PDFs, but now you need to install some really suspicious spyware-looking plugin for Adobe Reader to enforce some kind of DRM.

----

Anyway, hopefully I haven't distracted too much from the discussion at hand. I'm observing because I'm curious about collision handling too. One problem that is really bothering me is the simulation of Newton's Cradle.

Share on other sites
Quote:
Original post by Emergent
Quote:
 Original post by NumsgilGaussian Elimination, from my understanding, solves for the inverse of the matrix and then uses it to solve for the unknowns.

Thanks for your reply, Numsgil. But FYI -- and, apologies for this somewhat off-topic parenthetical discussion in this thread -- this is not correct. Solving a single system of equations with Gaussian elimination does not involve constructing an inverse. Just the opposite: You compute the inverse with GE by solving n equations (one for each element of the natural basis).

Ah, sorry. I haven't played with Gaussian Elimination since high school. Sort of skipped the section on that in Numerical Recipes and went straight to LU.

Quote:
 Anyway, hopefully I haven't distracted too much from the discussion at hand. I'm observing because I'm curious about collision handling too. One problem that is really bothering me is the simulation of Newton's Cradle.

It's what I'm playing with right now. I'm pretty sure you can construct a solver for newton's cradle using just local information, and solve the system globally at one time. Here's a powerpoint of what I've done so far (you might need texpoint to see the LaTeX render properly). It's incorrect in the slides, though, that as presented it will properly solve Newton's cradle. The result will be physically correct, but basically underestimate the restitution of the bodies. My guess is that you need to modify the right hand side (the 'b' terms in the ppt) somehow.

Anyway, since you're working on the same thing give it a look over and let me know any thoughts. That goes for anyone else interested as well. The 'j' terms in the ppt are the same (pretty sure) as the lambda terms for constraint solving. Just arrived at from a different direction. Likewise the Omega terms are basically J * M^-1 * JT, and the 'b' terms are basically J * v.