Please check logic behind LCP derivation

[XperienS]

Hello!
I'd be greatful if you'll check following statements about how to derive LCP formulation of rigid body dynamics and correct them if they're wrong.
The starting point of this derivation is D'Alembert-Lagrange variational principle, which states shortly that work of ideal constraint force in system is zero. So we must minimize constraints work to achieve ideal constraints:

min F(r) = F_c dot r
subject to c_e(x+r) = 0, c_i(x+r) >= 0;
where F_c - constraint force, r is displacement vector, c(r) is constraint function. I should mention here that in bilateral constraints non-zero work will cause constraint violation and in unilateral constraints valid displacement will not cause any constraint force. Also, later in this model could be introduced non-ideal constraints (e.g. friction), but there will be complementary constraints which restrict force and hence minimization to some finite number we need (-mu * N in case of Coulomb friction model).

To minimize this function, Lagrangian should be introduced:
L(r) = F(r) + lambda * (g - c(r)),
but we have g = 0, so
L(r) = F(r) - lambda * c(r);

Next we go straight to Karush-Kuhn-Tucker conditions:
Stationarity: dL / dr = dF(r) / dr - J^T * lambda = 0;
Primal feasibility: c_e(r) = 0, c_i(r) >= 0;

As we have dot product with constant in F(r), dF / dr = F_c, so we can rewrite first KKT condition:
F_c - J^T * lambda = 0;
but F_ext + F_c - M * a = 0, so
M * a - F_ext - J^T * lambda = 0;

Here we need to reformulate this in terms of lower-level derivatives (according to Cline, there will be problems with second derivatives), we assume first-order integrator (i will denote W == M^-1):
a = (v' - v) / dt;
=> v' = v + dt * W * F_ext + dt * W * J^T * lambda;

After these manipulations, we need to take a close look at constraint equations. Since we have constraints c(x + r) = c(x') and we don't know that resulting x' (assume r = v' * dt), we should approximate it via Taylor series in x:
c(x') = c(x) + dc(x)/dr * (x' - x) = c(x) + J * r = c(x) + dt * J * v';
We can substitute approximation and v' in primal feasibility (consider equality and inequality constraints mixed):
c(x) + dt * J * v' >= 0;
dt * J * W * J^T * lambda >= -c(x) / dt - J * (v + dt * W * F_ext);

Which is final LCP formulation.

[fishboy82]

Looks greate, and it's my first time to see using work minimization theory to deduce rigid body formulation,seems very short and elegant,at least I can't find any error in the deduce process , May be the only subtle problem is KKT need to promise constraint function to be continuously differentiable ,but this seems can be proved for every type of constraint.
Another question is how can suppose F_c is constant or irrelevant to diplacement r?

[XperienS]

fishboy82
Thanks for quick reply!

> Another question is how can suppose F_c is constant or irrelevant to diplacement r?
Well, this statement goes from that fact that when we minimizing work, we are choosing from virtual displacements, and virtual displacement is displacement which is determined by constraints at current time, hence constraint force is considered constant. You can see that in further substitution F_c + F_ext - M * a = 0, and in fact that a (not a') is approximated as (v' - v) / dt.

> it's my first time to see using work minimization theory to deduce rigid body formulation
Can you refer me to another type of derivation you meant?

[FD]

> Another question is how can suppose F_c is constant or irrelevant to diplacement r?

Well this is in fact a rather weak point, but it is obviously true if we look at the model that is used. We integrate the whole system in discrete steps, assuming that forces and velocities are constant during the timestep. Constraint forces calculated at some point in time are not changing untill the next simulation step; they are as discreet as everything else. This can be observed from the behaviour of dynamics engines, and I think this must be one of the assumptions.

What Xperiens says is that we actually use virtual displacements to advance the system. Virtual displacements are those that are allowed by constraints configuration at some given moment; the minimization functions actually choses the virtual displacement with minimal work of constraint forces, and then this (virtual) displacement is actually used to advance the system; since we consider virtual displacements, we only account for constrain forces corresponding to constraints configuration at the beginning of the frame, hence these forces are constant over the entire frame. for this to be true, we again must use the assumption that we advance the system by one of the virtual displacements, e.g. consider the contsraints being stationary at any moment. This is a more complex explaination of the same thing - we actually ASSUME these forces (as well as other forces and velocities) being contsnt over the time frame.

[fishboy82]

Thanks for FD and XperienS ,But I am still confused about F_c. As XperienS mentioned
F_c + F_ext - M * a = 0 and a = (v' - v) / dt and r = v' * dt
So we can easily get F_c + F_ext = M*(r/dt-v)/dt, But if we take a look at this formula
it seems F_c is obviously relevant to r, So is the assumption dF / dr = F_c right?

[XperienS]

fishboy82
Well, that is actually a VERY weak point of this kind of proof, but at first glance i can say in addition to what i've already said above that we're following this simple scheme for frame update:
1) Frame 0 state: x_0, v_0;
2) Calculate F_c from Frame 0 state: F_c(x_0, v_0);
3) Calculate a = W * F_ext + W * F_c; [equals to F_ext + F_c - M * a = 0];
4) Advance to Frame 1:
- v_1 = v_0 + a * dt;
- x_1 = x_0 + v_1 * dt = x_0 + r;
5) Frame 1 state: x_1, v_1;
It is obvious from this point that F_c doesn't actually relevant to r. All i wrote above is a bit tricky though

[XperienS]

To avoid this mess with F_c in explanation, i think we should go in more general way; and if we want to minimize work of constraint forces, we should take initial formulation, such as
F(r) = I(F_c * dr);
I denotes integral. Also, we have our discretization and linearization scheme, so we actually don't bother with that is happening between frames, so it seems that mentioned derivation is still correct:
dF(r)/dr = F_c;
I think that assumption of linear and discrete state propagation allows us to think that F_c is constant, but i'm not actually sure about that.

[fishboy82]

Thanks a lot ,I will take a more close look