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.

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.

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.

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.)