Решая LCP - Автор XperienS

Данная статья является обобщением знаний, полученных из работ по этой теме ( список работ представлени в конце статьи ); это выборка основных моментов, нужных для понимания основных принципов, почти без математики.

Решая LCP
LCP - Linear Complementarity Problem ( проблема линейного дополнения ) возникает тогда, когда нам нужно решить систему следующего вида:
b = Ax + q,
где b, x, q - векторы n x 1, A - матрица n x n;
действуют также так называемые дополнительные ограничения ( complementarity constraints ).
LCP нужно решать например если в симуляторе динамики твердого тела есть так называемые односторонние ограничения ( unilateral constraints ), такие как контактное соединение ( contact joint ) - сила должна действовать только тогда, когда тела нужно раздвинуть, но не должна их приближать. Для решения подобных систем существует два вида алгоритмов:
1. Прямые методы решения LCP;
Lemke: [ Michael Cline ]
Dantzig: [ Davig Baraff ]
2. Итеративные методы решения LCP;
Jacobi: [ Kenny Erleben ]
Gauss-Seidel: [ Erin Catto, Kenny Erleben ]

В этой заметке я остановлюсь лишь на итеративных методах так как за последнее время я удостоверился, что прямые методы лучше не использовать для игровых приложений. Они точнее, но и гораздо нестабильнее.

Итеративные методы решения LCP
Суть итеративных методов состоит в том, что ему дается некое начальное значение искомой переменной ( initial guess ), и после каждой итерации искомая переменная приближается от начального предположения к действительному решению. Очевидное преимущество - быстрота такого подхода по сравнению с прямыми методами решения, хотя и не такая высокая точность. Также можно ускорить работу, сделав так называемые тесты на сходимость ( convergence testing ), так же называемые критериями остановки ( stopping criteria ). Это позволит завершить работу солвера если требуемая точность достигнута. Этот способ будет описан позже.
Сначала я рассмотрю методы решения линейных систем вида Ax = b, без всяких дополнительных условий. Затем я покажу как один из методов решения линейных систем развить до метода решения MLCP ( Mixed LCP - присутствуют как ограниченные строки, так и неограниченные ). Итак, простейший итеративный метод решения линейных систем - Jacobi.
Математический вид данного метода:

x'_i = (b_i - S(L_i,j*x_j)_j=0,i-1 - S(U_i,j*x_j)_j=i+1,n-1) / A_i,i;

где S - сумма по переменной указанной мельче в пределах, указанных после равно; x' - переменная на текущем шаге; x - переменная на предыдущем шаге.
Код:
(
VecN - вектор с варьирующейся размерностью;
MatMN - матрица с изменяемой размерностью, причем m - кол-во строк, а n - кол-во столбцов.
)

VecN Solve_Jacobi(MatMN &A, VecN &b, int kMax = 5)
{
  if (A.m != b.n) return b;

  int i, j;
  int n = A.m;    //  Number of A's rows ( also size of vector b )

  VecN x = b;    // Initial Guess
  VecN y(n);
  float delta;

  // Jacobi Solver
  for (int k = 0; k < kMax; ++k)
  {
    for (i = 0; i < n; ++i)
    {
      delta = 0.0f;
      for (j = 0;   j < i; ++j) delta += A[i][j]*x[j];
      for (j = i+1; j < n; ++j) delta += A[i][j]*x[j];

      y[i] = (b[i] - delta) / A[i][i];
    }
    x = y;
  }

  return x;
}

Вышеуказанный метод является "параллельным"( это видно из-за того, что при вычислениях не используется x', а так же использование доп. вектора y: все переменные могут обновляться одновременно, потому что вычисления зависят только от результатов вычислений на предыдущем шаге ), в то время как метод Gauss-Seidel является "последовательным"( вычисления зависят как от результатов на текущем, так и на предыдущем шагах ):

x'_i = (b_i - S(L_i,j*x'_j)_j=0,i-1 - S(U_i,j*x_j)_j=i+1,n-1) / A_i,i;

Код, реализующий Gauss-Seidel solver:

VecN Solve_GaussSeidel(MatMN &A, VecN &b, int kMax = 5)
{
  if (A.m != b.n) return b;

  int i, j, n = A.m;

  VecN x = b;
  float delta;

  // Gauss-Seidel Solver
  for (int k = 0; k < kMax; ++k)
  {
    for (i = 0; i < n; ++i)
    {
      delta = 0.0f;

      for (j = 0;   j < i; ++j) delta += A[i][j]*x[j];
      for (j = i+1; j < n; ++j) delta += A[i][j]*x[j];

      x[i] = (b[i] - delta) / A[i][i];
    }
  }

  return x;
}

Несмотря на то, что метод Gauss-Seidel - последовательный, он более предпочтителен, нежели метод Jacobi. На практике он показывает более высокую скорость вычислений и стабильность.

SOR
Так же для увеличения стабильности данного метода используется так называемая последовательная верхняя релаксация ( successive overrelaxation - SOR ). Для того, чтобы понять, как она работает, нужно определить остаточный вектор( residual vector ) r:
остаточный вектор для x в данной системе: r = Ax - b;
тогда получим, что
A_i,i*x_i + r'_i,i = A_i,i*x'_i;
переставив местами:
x'_i = x_i + r'_i,i/A_i,i;
для более быстрой сходимости метода вводят новую переменную w ( некоторые ее значения будут обеспечивать более быструю сходимость ), называемую параметром релаксации ( в данном случае для метода Gauss-Seidel ):
x'_i = x_i + w * r'_i,i/A_i,i;
теперь метод преобразовался в SOR Gauss-Seidel, параметр w должен быть 1 + 0.2 ( когда 0 < w < 1 это нижняя релаксация - under-relaxation, а w > 1 это верхняя релаксация; у меня в данный момент w = 0.9 ) и код для решения линейной системы методом SOR-GS:

VecN Solve_SOR(MatMN &A, VecN &b, float w = 0.9f, int kMax = 5)
{
  if (A.m != b.n) return b;

  int i, j, n = A.m;

  VecN x = b;
  float delta;

  // Gauss-Seidel with Successive OverRelaxation Solver
  for (int k = 0; k < kMax; ++k)
  {
    for (i = 0; i < n; ++i)
    {
      delta = 0.0f;

      for (j = 0;   j < i; ++j) delta += A[i][j]*x[j];
      for (j = i+1; j < n; ++j) delta += A[i][j]*x[j];

      delta = (b[i] - delta) / A[i][i];
      x[i] += w * (delta - x[i]);
    }
  }

  return x;
}

Тест на сходимость или Критерий остановки
Вообще-то по-хорошему надо, чтобы итеративные метода продолжали решать, пока тест на сходимость не пройден. Но в действительности этот тест может быть использован в купе с максимальным количеством итераций, чтобы ограничить количество итераций ( а следовательно, ускорить работу алгоритма ). Для того, чтобы провести тест на сходимость, нужно ввести некую нужную точность CONV_EPS, тогда тест на сходимость будет выглядеть примерно так:
пусть x' - решение на текущем шаге, x - решение на предыдущем и определена точность CONV_EPS > 0
||x' - x|| < CONV_EPS;
где || || - векторная норма.
Один из способов проверки на сходимость это просто проверка величины delta. Если delta < CONV_EPS для всех i на данной итерации, то тест на сходимость пройден и можно завершать решение системы. У меня пока этот тест не реализован, но сделать его, особенно по второму методу не составит особого труда.

Страницы: 1 2 Следующая »

#dynamics, #LCP, #физика, #rigid body

15 апреля 2006 (Обновление: 21 янв 2010)

Комментарии [2]