Linear Complementarity Problem (LCP)

The core of an LCP is the complementarity condition, expressed as 0 ≤ z ⟂ w ≥ 0. This means the solution variables z and w must be non-negative and orthogonal: for each component i, either z_i = 0 or w_i = 0 (or both). In physics, this elegantly encodes the non-penetration and non-adhesion of contact: a contact force (z) can only be positive if the bodies are exactly touching (w=0), and if bodies are separating (w>0), the contact force must be zero.

The 'Linear' in LCP refers to the relationship w = Mz + q, where M is a known matrix and q is a known vector. This linear system is derived from the physics of the problem. For contact, M encodes the system's mass matrix and contact geometry (Jacobians), while q contains external forces and relative velocities. This structure allows the problem to be solved with efficient, specialized numerical algorithms rather than general non-linear solvers.

LCPs naturally model unilateral constraints, which are inequalities that can be active or inactive. A classic example is contact: bodies cannot push into each other (the normal force must be ≥ 0), but they can freely separate. This is distinct from bilateral constraints (like a hinge joint), which are always active equalities. The LCP's solution automatically determines which contacts are active (exerting force) and which are inactive (separating), a process known as contact state determination.

Coulomb friction is non-linear, but it can be incorporated into an LCP through a polyhedral approximation of the friction cone. The continuous cone of allowable friction forces is approximated by a pyramid with 4, 8, or more sides. Each side becomes a linear complementarity constraint. This transforms the non-linear friction problem into a larger, but still linear, complementarity problem that can be solved with the same numerical machinery, enabling stable simulation of sticking and sliding.

A standard algorithm for solving LCPs is Lemke's algorithm, a pivotal method similar to the simplex algorithm for linear programming. It systematically explores vertices of a complementary cone to find a solution. For physics simulation, the Projected Gauss-Seidel (PGS) method is often preferred for its simplicity and efficiency in handling large, sparse systems, even though it provides an approximate solution. The choice of solver balances accuracy, speed, and robustness for real-time applications.

The LCP is intimately related to convex optimization. Solving the LCP 0 ≤ z ⟂ Mz + q ≥ 0 is equivalent to finding a Karush-Kuhn-Tucker (KKT) point for a specific Quadratic Program (QP) when M is positive semi-definite. This link provides a rich theoretical foundation and allows the use of optimization-based techniques like Sequential Quadratic Programming (SQP) for more complex, non-linear contact problems, bridging the gap between game physics and high-precision engineering simulation.

A constraint solver is the algorithmic core of a physics engine that calculates the forces or impulses necessary to satisfy a set of constraints—such as contact non-penetration, joint limits, or friction—within a numerical tolerance. The LCP is a specific mathematical framework often used to formulate and solve these constraints. The solver iteratively adjusts velocities and positions to prevent interpenetration while respecting physical laws like Coulomb friction.

• Primary Role: Enforces physical rules (no sinking, sliding limits) in a simulated world.
• Iterative Nature: Typically uses algorithms like Projected Gauss-Seidel (PGS) to find a feasible solution.
• Connection to LCP: Many modern solvers frame contact and friction constraints as an LCP or a related optimization problem to find a globally consistent set of forces.

Coulomb friction is the classical, dry friction model used in physics simulation. It defines two regimes:

• Static Friction: The force that resists the initiation of sliding, with a maximum magnitude proportional to the normal force.
• Kinetic Friction: A constant opposing force during sliding, also proportional to the normal force.

The LCP is instrumental in modeling Coulomb friction because it naturally encodes the complementarity conditions between force and motion:

• If there is no sliding (relative tangential velocity is zero), the friction force can be anywhere within its allowable friction cone.
• If sliding occurs, the friction force must be at the cone's boundary, opposing the direction of slip. The LCP solver finds the exact force that satisfies these either-or conditions at every contact point.

A friction cone is the geometric representation of all allowable friction forces at a contact point, as defined by Coulomb's law. It is a cone whose:

• Axis is aligned with the contact normal.
• Angle is given by arctan(μ), where μ is the coefficient of friction. Any resultant tangential force must lie within or on the surface of this cone to be physically admissible.

In LCP formulations, the friction cone is often linearized into a polyhedral approximation (e.g., a pyramid with 4 or 8 edges) to transform the nonlinear friction constraints into a set of linear complementarity conditions. This allows the use of efficient, dedicated LCP solvers to determine if a contact is sticking (force inside cone) or sliding (force on cone's edge).

The penalty method is an alternative to LCP-based constraint resolution. Instead of solving for exact non-penetration constraints, it models contacts as stiff spring-damper systems.

• Mechanism: When two bodies penetrate, a spring force proportional to the penetration depth pushes them apart. A damping term proportional to penetration velocity adds stability.
• Contrast with LCP: The LCP is a hard constraint method that strictly enforces no penetration at the velocity level. The penalty method is a soft constraint that allows small, numerical penetrations to generate a counteracting force.
• Use Cases: Penalty methods are simpler to implement and can be more stable for very stiff, non-interacting constraints, but they can cause "bouncing" and are less accurate for complex, simultaneous contacts compared to LCP solvers.

Impulse-based dynamics is a simulation method that resolves collisions and contacts by applying instantaneous changes in velocity (impulses) rather than computing continuous forces over a time step.

• Process: When a collision is detected, an impulse is calculated and applied to instantly change the bodies' velocities to satisfy conservation laws and collision properties (like restitution).
• Relation to LCP: Sequential impulse-based solvers, a popular variant, are mathematically similar to solving an LCP using the Projected Gauss-Seidel method. Each iteration applies impulses to satisfy one constraint at a time, progressively converging to a solution that respects all contacts and friction conditions simultaneously.
• Advantage: Known for stability and efficiency in real-time applications like video games, where exact physical accuracy is less critical than robust performance.

Baumgarte stabilization is a constraint stabilization technique used to control numerical drift in simulations that use constraint-based methods like the LCP. Constraint drift is the accumulation of small numerical errors causing bodies to slowly penetrate or joints to separate.

• Function: It modifies the constraint equation by adding a corrective feedback term. For a position constraint C(p)=0, the solver is told to achieve Ċ + βC = 0, where β is a damping parameter. This term acts like a virtual spring pulling the system back to a valid state.
• Role in LCP: While the LCP is often formulated at the velocity level, Baumgarte stabilization provides a practical way to incorporate position-level correction, ensuring long-term simulation stability without requiring perfectly solved, drift-free constraints at every infinitesimal step.