XPBD: розширена позиційна динаміка

1. Чому звичайного PBD було недостатньо

Класична Position-Based Dynamics (Müller та ін., 2007) симулює обмеження — відстань, згин, об'єм — шляхом прямого проєктування позицій частинок так, щоб функція обмеження C(x) = 0 виконувалася, а потім повторює це проєктування фіксовану кількість разів за кадр (у стилі Gauss-Seidel). Це швидко, безумовно стабільно і тривіально реалізується. Але є проблема: у самому кроці проєктування немає нічого, що відповідало б реальному значенню жорсткості в Н/м. "Жорсткість", яку ви отримуєте, — це артефакт двох непов'язаних чисел: кількості ітерацій розв'язувача та кроку часу, помножених одне на одне.

Зменшіть крок часу вдвічі або подвойте кількість ітерацій — і мотузка, що раніше помітно розтягувалась, раптом починає поводитися як жорсткий стрижень. Це робило PBD непридатним для всього, що мало відповідати реальній поведінці матеріалу, і означало, що налаштована на одній машині симуляція тканини могла виглядати абсолютно інакше на іншій машині з іншою частотою кадрів.

2. Комплаєнс і множники Лагранжа

XPBD переформульовує проєктування як розв'язання задачі умовної оптимізації прямо з лагранжевої механіки. Кожне обмеження C_j(x) = 0 отримує значення комплаєнсу α_j — обернену величину фізичної жорсткості, виміряну в м/Н (метрах на Ньютон, тобто наскільки обмеження "поступається" на одиницю сили):

Повна потенційна енергія обмеження дорівнює U(x) = ½ · α⁻¹ · C(x)². Мінімізація цієї енергії за умови рівнянь руху вводить множник Лагранжа λ для кожного обмеження — фізично λ є накопиченим імпульсом сили обмеження. Ключовий трюк XPBD полягає в тому, щоб оновлювати λ інкрементально разом з оновленням позиції, замість розв'язання повної лінійної системи:

Задайте α = 0 — і це точно зводиться до формули проєктування класичного PBD: XPBD є строгою надмножиною. Задайте α > 0 — і ви отримаєте пружину з фізично осмисленою жорсткістю, незалежною від кількості ітерацій.

3. Цикл розв'язувача XPBD

Повний алгоритм на кадр чергує субкроки, кожен з яких прогнозує позиції, розв'язує всі обмеження один раз (з обнуленням λ на початку субкроку, а не кадру) та виводить швидкість зі зміни позиції:

function xpbdStep(particles, constraints, dt, substeps) {
  const h = dt / substeps; // крок субкроку

  for (let s = 0; s < substeps; s++) {
    // 1. Прогноз позицій із зовнішніми силами (гравітація, вітер...)
    for (const p of particles) {
      if (p.invMass === 0) continue; // закріплена/статична
      p.vel.y += -9.8 * h;
      p.prevPos.x = p.pos.x; p.prevPos.y = p.pos.y;
      p.pos.x += p.vel.x * h;
      p.pos.y += p.vel.y * h;
    }

    // 2. Скидання множників Лагранжа для цього субкроку
    for (const c of constraints) c.lambda = 0;

    // 3. Розв'язати кожне обмеження один раз (прохід Gauss-Seidel)
    for (const c of constraints) solveConstraint(c, h);

    // 4. Вивести швидкість зі зміни позиції (основна ідея PBD)
    for (const p of particles) {
      if (p.invMass === 0) continue;
      p.vel.x = (p.pos.x - p.prevPos.x) / h;
      p.vel.y = (p.pos.y - p.prevPos.y) / h;
    }
  }
}

Зверніть увагу — тут немає окремого проходу "розв'язання швидкості", як у рушіях на основі імпульсів: швидкість завжди є побічним продуктом зміни позиції, поділеної на h. Саме це робить PBD і XPBD безумовно стабільними: ви ніколи не можете розійтися до нескінченної швидкості, бо швидкість обмежена тим, наскільки далеко позиція реально може змінитись за один субкрок.

4. Виведення обмеження відстані

Робоче обмеження для тканин і мотузок — це обмеження відстані, яке утримує дві частинки на відстані спокою L₀:

function solveDistanceConstraint(c, h) {
  const { p1, p2, restLength, compliance } = c;
  const dx = p1.pos.x - p2.pos.x;
  const dy = p1.pos.y - p2.pos.y;
  const dist = Math.hypot(dx, dy);
  if (dist < 1e-9) return;

  const C = dist - restLength;
  const nx = dx / dist, ny = dy / dist;

  const alphaTilde = compliance / (h * h);
  const denom = p1.invMass + p2.invMass + alphaTilde;
  if (denom < 1e-9) return;

  const dLambda = (-C - alphaTilde * c.lambda) / denom;
  c.lambda += dLambda;

  p1.pos.x +=  p1.invMass * dLambda * nx;
  p1.pos.y +=  p1.invMass * dLambda * ny;
  p2.pos.x -=  p2.invMass * dLambda * nx;
  p2.pos.y -=  p2.invMass * dLambda * ny;
}

Той самий скелет — обчислити C, обчислити ∇C для кожної частинки, розв'язати відносно Δλ, застосувати зважені зміщення позицій — працює і для обмежень згину (двогранний кут між двома трикутниками), обмежень об'єму (знаковий об'єм тетраедра) та обмежень площі. Змінюється лише C та її градієнт.

5. Чому субкроки, а не ітерації розв'язувача

Оригінальна стаття PBD фіксувала крок часу на рівні частоти кадрів (наприклад, 1/60 с) і збільшувала кількість ітерацій для покращення збіжності та жорсткості. Автори XPBD виявили, що це погано збігається — прохід Gauss-Seidel по багатьох обмеженнях поширює корекцію лише на одну ланку за ітерацію, тож довгій смузі тканини потрібні десятки ітерацій, щоб відчуватись натягнутою.

Натомість XPBD надає перевагу багатьом малим субкрокам з одним розв'язанням обмежень кожен, замість одного великого кадру з багатьма ітераціями. За той самий обчислювальний бюджет субкрокінг збігається значно швидше, бо член комплаєнсу α̃ = α / h² зменшується разом з h — менші субкроки автоматично роблять кожне обмеження жорсткішим і точнішим, без переналаштування жодного параметра.

6. Демпфування без багів втрати енергії

Наївне демпфування швидкості після розв'язання (v *= 0.98) видаляє енергію нерівномірно і може знову зробити симуляцію залежною від частоти кадрів. Натомість XPBD рекомендує демпфування на основі обмежень: додати невелику додаткову корекцію позиції, пропорційну відносній швидкості вздовж градієнта обмеження, масштабовану комплаєнсом демпфування β та субкроком h:

Це утримує демпфування пропорційним фактичній відносній швидкості на обмеженні (як фізичний демпфер/амортизатор), тож симуляція тканини, запущена на 20 субкроках, і та, що на 30, осідають до однакової форми спокою з однаковою швидкістю.

7. Повна реалізація: підвісний ланцюг

Зберемо все разом — ланцюг частинок, з'єднаних обмеженнями відстані, з першою частинкою закріпленою (invMass = 0):

class Particle {
  constructor(x, y, invMass = 1) {
    this.pos = { x, y };
    this.prevPos = { x, y };
    this.vel = { x: 0, y: 0 };
    this.invMass = invMass;
  }
}

class DistanceConstraint {
  constructor(p1, p2, compliance = 0) {
    this.p1 = p1; this.p2 = p2;
    this.restLength = Math.hypot(p1.pos.x - p2.pos.x, p1.pos.y - p2.pos.y);
    this.compliance = compliance; // 0 = жорстке, більше = розтяжніше
    this.lambda = 0;
  }
}

// Побудова ланцюга з 12 ланок, перша ланка закріплена на місці
const chain = [];
for (let i = 0; i < 12; i++)
  chain.push(new Particle(i * 0.2, 2, i === 0 ? 0 : 1));

const links = [];
for (let i = 0; i < chain.length - 1; i++)
  links.push(new DistanceConstraint(chain[i], chain[i + 1], 0.0001)); // злегка податлива сталь

// Головний цикл (викликати раз на відрендерений кадр, dt ≈ 1/60)
function simulate(dt) {
  xpbdStep(chain, links, dt, 24); // 24 субкроки
}

8. Застосування: тканини, стержні та м'які тіла

• Тканина: сітка обмежень відстані (структурні + зсувні + згинальні) з комплаєнсом для кожного обмеження дає реалістичне драпірування без артефакту "PBD-повітряної кульки" надмірної жорсткості, характерного для класичних реалізацій.
• Стержні та волосся: XPBD природно поширюється на обмеження, засновані на орієнтації (з використанням кватерніонів як узагальнених координат), для жорсткості згину/скручування в моделях стержнів Коссера — саме так Blender симулює пасма волосся.
• М'які тіла: тетраедральні сітки з обмеженнями збереження об'єму та обмеженнями енергії деформації типу neo-Hookean, розв'язані тим самим оновленням Δλ, відтворюють FEM-подібну пружну поведінку за частку вартості.
• З'єднання жорстких тіл: формулювання через комплаєнс узагальнюється на шарнірні з'єднання з 6 ступенями свободи (шарнір, кульовий шарнір, призматичне з'єднання), тому і PhysX 5, і фізичний стек Unity перевели свої розв'язувачі з'єднань на оновлення у стилі XPBD.

9. Типові пастки

• Скидання λ у неправильній області видимості: λ має скидатися до нуля раз на субкрок, а не раз на кадр і не раз на ітерацію Gauss-Seidel всередині субкроку — помилка тут непомітно ламає всю математику комплаєнсу.
• Забути про α̃ = α / h²: використання сирого значення комплаєнсу замість перерахованого на субкрок робить жорсткість знову залежною від кількості субкроків, зводячи нанівець весь сенс XPBD.
• Нульовий комплаєнс без захисту від нульового знаменника: коли обидві частинки в обмеженні закріплені (invMass = 0 для обох), знаменник дорівнює нулю — завжди перевіряйте ділення на нуль перед обчисленням Δλ.
• Артефакти, залежні від порядку: проходи Gauss-Seidel послідовні, тож порядок розв'язання обмежень трохи впливає на результат; для паралельних (GPU) розв'язувачів використовуйте розфарбування графа, щоб згрупувати неконфліктуючі обмеження та розв'язувати їх одночасно без гонок даних.