Constraint satisfaction у фізиці: метод множників Лагранжа

1. Чому не просто пружини?

Найпростіший спосіб «прив'язати» дві точки — з'єднати їх дуже жорсткою пружиною: F = −k·(|r| − L)·r̂. Проблема в тому, що жорсткість k, потрібна для «нерозтяжної» мотузки, робить систему чисельно жорсткою (stiff) — явний інтегратор вимагає надзвичайно малого dt, інакше симуляція вибухає осциляціями.

Метод множників Лагранжа дозволяє знайти точну силу, яка утримує обмеження — без будь-якої жорсткості чи власної частоти коливань. Це та сама математика, яка лежить в основі узагальненої механіки Лагранжа: мінімізація дії за умови обмежень.

2. Механіка Лагранжа з обмеженнями

Обмеження записують як неявну функцію C(q) = 0, де q — узагальнені координати системи. Наприклад, для жорсткого стрижня довжини L між точками p₁ і p₂:

Рівняння руху з обмеженнями (constrained Euler-Lagrange) додають до звичайних сил доданок λ · ∇C — силу, спрямовану вздовж градієнта обмеження, з невідомим масштабом λ (множник Лагранжа):

Множник λ обчислюється так, щоб обмеження точно виконувалось на кожному кроці — не наближено, як з пружиною, а точно (з точністю до похибки чисельного інтегрування).

3. Що таке множник λ насправді

λ — це величина сили обмеження, необхідної, щоб утримати систему на многовиді C(q) = 0. Фізично це натяг мотузки, реакція опори, сила зчеплення в шарнірі. Знак λ несе фізичний зміст: для мотузки λ ≥ 0 (мотузка тягне, але не штовхає) — це однобічне обмеження (unilateral constraint), на відміну від жорсткого стрижня, де λ може бути будь-якого знаку (двобічне, bilateral).

4. Обмеження на рівні швидкості

У фізичних рушіях реального часу зручніше працювати не з позиційним обмеженням C(q) = 0, а з його похідною за часом — обмеженням на швидкість:

Підставляючи рівняння руху, отримуємо лінійну систему відносно λ:

Матрицю J · M⁻¹ · Jᵀ називають ефективною масою (effective mass) обмеження — вона показує, наскільки «легко» система реагує на одиничний імпульс уздовж напрямку обмеження.

5. Baumgarte-стабілізація

Робота лише з Ċ = 0 утримує обмеження на рівні швидкості, але не виправляє накопичену позиційну похибку (drift) через похибки чисельного інтегрування — мотузка поступово «розтягується» на десяті частки відсотка за кроками. Baumgarte-стабілізація додає корекційний доданок, пропорційний позиційній похибці:

Параметр β ∈ [0.1, 0.3] контролює швидкість корекції дрейфу: занадто велике β «висмикує» систему занадто різко (додає штучну енергію), занадто мале — дрейф накопичується видимо. Це компроміс між точністю і стабільністю — типова інженерна проблема числових методів.

6. Sequential impulse solvers

Коли обмежень багато (десятки контактів, шарнірів у ragdoll), розв'язання повної системи J·M⁻¹·Jᵀ·λ = b напряму — дорого. Практичні рушії (Box2D, Bullet, PhysX) використовують Sequential Impulse (Projected Gauss-Seidel): розв'язують кожне обмеження окремо, ітеративно, приймаючи інші за фіксовані, поки система не збіжиться.

Зазвичай достатньо 4–10 ітерацій для правдоподібної збіжності в реальному часі — точний розв'язок не потрібен, важлива стабільність кадр за кадром.

7. Порівняння методів обмежень

8. Псевдокод

Обмеження відстані з Baumgarte (JS)

// Розв'язати λ для обмеження відстані з Baumgarte-корекцією
function solveDistanceConstraint(a, b, restLength, dt, beta = 0.2) {
  const delta = sub(b.pos, a.pos);
  const dist = length(delta);
  const n = scale(delta, 1 / dist);        // одиничний напрямок
  const C = dist − restLength;               // позиційна похибка

  // Відносна швидкість вздовж напрямку обмеження
  const relVel = dot(sub(b.vel, a.vel), n);

  // Ефективна маса вздовж n
  const invMassSum = a.invMass + b.invMass;
  const bias = (beta / dt) * C;             // Baumgarte-доданок

  const lambda = −(relVel + bias) / invMassSum;

  // Застосувати імпульс Jᵀ·λ до обох тіл
  a.vel = sub(a.vel, scale(n, lambda * a.invMass));
  b.vel = add(b.vel, scale(n, lambda * b.invMass));
}