Введение

Представлен метод формирования эффективных в вычислительном плане уравнений динамики манипуляционных роботов.

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

В разделе 2 приводится формализация метода последовательного формирования систем координат звеньев манипулятора для описания его кинематики. Для расчета кинематических и динамических величин используются матрицы преобразования координат размера 3х3 и век­тора относительных перемещений. Показано применение метода для бортового манипулятора космического корабля “Буран” и промышленного робота-манипулятора РМ-01. В разделах 3 – 5, представлен метод формирования уравнений динамики манипуляторов в форме уравнений Лагранжа II рода. Метод применим для манипуляторов с вращательными и поступательными шарнирами, соседние оси которых перпендикулярны или параллельны. Приведена оценка вычислительных затрат, необходимых для расчета кинематики и динамики манипуляторов данного типа.

В разделе 6 представлено использование языка символьных вычислений REDUCE для сим­вольного вывода уравнений динамики манипуляторов. Приводится анализ вычислительной эффективности уравнений в символьном виде для манипуляторов с 2 и 3 степенями свободы.

Для замыкания системы уравнений динамики необходимо получить выражения для обобщенных моментов в шарнирах. Они определяются типами и параметрами двигателей, механических передач (редукто­ров), а также особенностями системы управления робота. В разделе 7 рассмотрены электромеханические привода общего вида с двигателями пос­тоянного тока, обратимыми зубчатыми редукторами и замкнутыми по положению и скорости следящими системами. Представлены различные по сложности модели приводов, в которых учитываются упругость в шарнирах и нелинейные элементы: люфт, сухое трение, муфты пре­дельного момента, тормоза.

1. Обзор методов описания кинематики и динамики манипуляционных роботов

Формирование эффективных уравнений динамики манипуляционных роботов, которые могут быть рассчитаны на ЭВМ за минимальное время, является одной из важнейших задач в робототехнике. Ее решение необходимо для моделирования динамики манипуляторов в масштабе реального времени, для разработки эффективных алгорит­мов управления роботами с учетом динамики [1], для повышения эффективности исследования и разработки манипуляторов.

Одни из первых результатов в этой области принадлежат Кейну [2] и Виттенбургу [3]. Полученные ими уравнения справедливы не только для роботов, но и для более широкого класса систем, состоящих из шарнирно связанных твердых тел. В дальнейшем было раз­работано большое количество алгоритмов формирования динамических уравнений манипуляторов, в которых использовались различные способы описания кинематики, расчета кинематических и динамических величин, а также различные формы уравнений динамики системы тел.

Описание кинематики – это способ задания систем координат, связанных со звеньями манипулятора, и выбора параметров, которые однозначно определяют взаимное положение звеньев и конфигурацию всего манипулятора. В представлении Денавита-Хартенберга [4] начала систем координат расположены в шарнирах, а их оси формируются по правилам, которые определяются кинематикой манипулятора. В другом методе описания кинематики [5-7] локальные сис­темы координат привязаны к центрам масс звеньев, а их оси направлены вдоль главных осей инерции. Параметры, определяемые относи­тельно таких систем координат, удобны для динамического анализа. В настоящей работе используется метод последовательного формиро­вания систем координат, предложенный в [8] (его описание при­водится в разделе 2).

Еще одной характеристикой методов математического моделирования манипуляторов является способ расчета кинематических и динамических величин, определяющих математическую модель манипу­лятора. Для этого используются однородные координаты и матрицы преобразования координат размерности 4x4, определяющие относи­тельное положение и ориентацию звеньев манипулятора [9-11]; матрицы поворотов размерности 3x3 и вектора относительных пере­мещений [9, 12, 13]; формулы Родриго (впервые применены в [7], далее использовались в [5, 14, 15]); ортогональные тензоры [16]; кватернионы [17]; метод векторных параметров с использованием групп Ли [18, 19].

Хотя вычислительная эффективность того или иного метода формирования динамических уравнений зависит в первую очередь от особенностей его реализации (использования рекурсивных преобразований, динамических аналогий и др.), можно отметить и сущест­венную роль выбора подходящего способа расчета модели манипуля­тора. Например, матрицы преобразования однородных координат раз­мерности 4x4, обладающие универсальностью в кинематическом опи­сании, практически не используются в задачах реального времени из-за больших вычислительных затрат, необходимых для выполнения операций над ними [1]. В то же время, использование матриц по­воротов размера 3x3 позволяет получить эффективные алгоритмы расчета кинематики и динамики, что показано в [1, 9] и в данной работе. Эффективно использование кватернионов, ортогональных тензоров (с их помощью получен самый быстрый алгоритм решения обратной задачи динамики [16]), однако в ряде задач (например, при управлении в декартовых осях, [20]) предпочтительнее исполь­зовать матричные представления.

При выводе уравнений динамики манипуляторов используются различные законы и формулировки общих уравнений динамики систем. Среди них можно выделить методы, основанные на уравнениях Лаг­ранжа, Ньютона-Эйлера, Д'Аламбера, Гаусса, Аппеля, Кейна.

Уравнения динамики в форме Лагранжа впервые были получены в работе Uicker [11] и получили дальнейшее развитие в плане повы­шения эффективности в работах Kahn [10], Vukobratovic [21], Mahil [14], Renaud [22], Thomas и Tesar [13]. Все перечисленные методы позволяли решать прямую и обратную задачу динамики, были удобны в алгоритмической реализации (кроме Renaud), но обладали низкой вычислительной эффективностью. Waters и Hollerbach [9] применили рекурсивные преобразования при выводе динамических уравнений, причем в [9] при использовании матриц поворотов 3x3 было получено значительное сокращение числа операций, но эти методы позволяли решать лишь обратную задачу динамики, поэтому не были пригодны для моделирования. Рекурсивные преобразования и формулы Родриго использовали Vukobratovic и Potkonjak [21], при­чем их метод позволял решать и прямую задачу динамики, хотя его вычислительная эффективность и не столь высока. Значительный прогресс в сокращении числа операций достигнут в работах Renaud [22] и Li [23], также применивших рекурсивные соотношения.

Среди методов, использующих уравнения Ньютона-Эйлера и позволяющих решать прямую задачу динамики, отметим работы Vukobra­tovic и Stepanenko [7], Orin и Walker [24], Armstrong [25], Wang и Ravani [15]. Во всех работах применяется рекурсия, причем в [7] в результате получаются нерекурсивные выражения для динамических коэффициентов, которые можно использовать для анализа динамики. Отметим также алгоритм Balafoutis [16], в котором за минимальное число операций решена обратная задача динамики.

Использование принципа Д'Аламбера для уравнений Лагранжа позволяет получить достаточно эффективные динамические соотношения, в которых в явном виде отражены эффекты влияния вращатель­ного и поступательного движения звеньев на динамику манипулятора [1, 26]. В работе Попова [27] получены уравнения динамики в яв­ном виде с использованием уравнений Аппеля, позволяющие решать прямую и обратную задачи динамики. Уравнения Кейна особенно эффективны для расчета обобщенных моментов манипуляторов с замкнутыми кинематическими цепями [28, 29]. В работах Коноплева [30] для описания динамики манипуляционных систем применяются агрегатив­ные модели; метод удобен для применения символьных преобразова­ний. Погореловым [31] разработан пакет программ моделирования динамики широкого класса механических систем, включая роботы-манипуляторы. Интересный подход и язык программирования уравнений движения сложных механических систем, состоящих их твердых тел, предложен Сазоновым [32].

Среди самых современных методов моделирования динамики манипуляторов отметим подходы, основанные на использовании нейронных  сетей [33, 34], пространственных операторов [35, 36], групп Ли [19], методов нечеткой логики [37]. Для описания динамики сложных структур (параллельных роботов, манипуляторов с большой избыточностью степеней подвижности, роботов-гуманоидов) используются методы расчета динамики в операционном пространстве роботов и другие. Полный анализ основных достижений в области моделирования динамики роботов, начиная с первых работ 60-70 годов прошлого века по 2000 г., дан Featherstone и Orin в [38].

В таблице 1 представлены результаты анализа методов описа­ния динамики      манипуляторов    по форме   уравнений,    вычислительной

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

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

В последние годы для повышения эффективности уравнений динамики широко применяются символьные преобразования [5, 39, 40, 41] и алгоритмы распараллеливания вычислений [42, 43, 44]. Применение символьных преобразований, как показано в [40], способствует уравниванию различных алгоритмов по вычисли­тельной эффективности. Поэтому важными критериями оценки алго­ритма становятся  хорошая  алгоритмизуемость (удобство программи­рования), замкнутость уравнений, возможность применения символь­ных преобразований и алгоритмов распараллеливания.

Кроме того, можно отметить, что в работах ряда авторов [9, 45] указывается, что во многих случаях (например в задачах управления роботами) наиболее подходящим способом описания дина­мики являются уравнения Лагранжа. Это отмечено и в работе [46], где утверждается, что при реализации динамических алгоритмов на параллельных процессорах зависимость данных в методах, основанных на уравнениях Ньютона-Эйлера, намного сильнее, чем при использовании уравнений Лагранжа.

Таблица 1.

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

2. Вывод основных кинематических соотношений

Рассмотрим n-звенный манипулятор с вращательными шарнирами. Для описания его кинематики необходимо с каждым из звеньев свя­зать систему координат, относительно которой будут рассчитывать­ся параметры звена, а также ввести обобщенные координаты, т.е. задать направления отсчета углов поворота звеньев.

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

С основанием манипулятора Т₀ свяжем неподвижную систему координат (СК) S₀, потребовав, чтобы одна из ее осей совпадала с осью первого поворота e_{rot 1}, а начало находилось в центре перво­го шарнира - точке О₁ (см. рис. 1).

Рисунок 1.                  Рисунок 2а.                        Рисунок 2б.

Теперь построим систему координат S₁, связанную с первым звеном Т₁. Ее начало выберем в центре первого шарнира. Одну из координатных осей СК S₁ направим по оси первого поворота. Если ось второго поворота e_{rot 2} || e_{rot 1} , то направления 2-х оставших­ся координатных осей выбираются с учетом геометрии звена (как правило, одна ось расположена вдоль звена, а другая в перпенди­кулярном направлении, см. рис. 2а). Если e_{rot 2} ^ e_{rot 1} , то одну из осей в плоскости, перпендикулярной e_{rot 1}, направим параллель­но вектору e_{rot 2}, а третья координатная ось будет дополнять выб­ранные две до правой ортогональной тройки (рис. 2б).

Таким образом, как СК S₀ и S₁, так и СК S₁ и S₂ будут иметь одну координатную ось с общим направлением, определяемым для S₀ и S₁  вектором e_{rot 1}, а для S₁  и S₂  вектором e_{rot 2}. Поэтому  матрицы перехода от S₀ к S₁ и от S₁ к S₂ будут матрицами  поворотов относительно соответствующих координатных осей : C_i Î {C_x, C_y, C_z}, где C_x, C_y, C_z - матрицы поворотов относительно осей x, y и z со­ответственно.

Аналогичным образом определим системы координат S₂, S₃, .., S_n оставшихся звеньев манипулятора (S_i связана со звеном T_i, а ее начало O_i находится в центре i-ого шарнира).

Для однозначного определения СК S₁, .., S_n необходимо задать направления отсчета обобщенных координат q₁, .., q_n (q_i определяет переход от S_i-1  к S_i). С этой целью введем вектора  e⁰_qi Î T_i-1 (направлен  по одной из координатных осей СК S_i-1, перпендикуляр­ных e_{rot i}) и  e¹_qi Î T_i (направлен по одной из координатных осей СК S_i, перпендикулярных e_{rot i}).

Рисунок 3.

Угол q_i будем отсчитывать от e_qi⁰ к e_qi¹; положительное направление отсчета определяется поворотом от e_qi⁰ к e_qi¹ относитель­но вектора e_{rot  i} против часовой стрелки.

Итак, построены системы координат звеньев и заданы обобщен­ные координаты манипулятора. Для полного кинематического описа­ния манипулятора осталось определить матрицы поворотов C_i, за­дающие ориентацию соседних систем координат S_i-1 и S_i, и вектора переноса l_i-1 = O_i-1O_i, определяющие сдвиг между ними (l₀ = 0, т.к.   O₀= O₁; l_n = O_nG, где G - конечная точка (схват) манипулятора).

С этой целью рассмотрим манипулятор в “начальном” состоя­нии, т.е. при q₁ = q₂ = ... = q_n = 0. Задав СК S₀ (т.е. оси координат х₀ , y₀ , z₀), выполним последовательно параллельные переносы S₀ в точки O₁ , O₂, ..., O_n. При этом будут определены оси систем коорди­нат S₁, ..., S_n (см. рис. 4в). Теперь легко для каждого звена опре­делить матрицы C_i и вектора l_i.

На рисунках 4а-4в иллюстрируется применение метода последовательного формирования систем координат к описанию кинематики манипулятора большого космического манипулятора (БКМ) – бортового манипулятора космического корабля “Буран”, рис. 4а.

Сперва рассмотрим манипулятор в общем положении и введем локальные СК и направления отсчета углов поворота q₁, q₂, ..., q_n, руководствуясь описанными выше правилами, а также используя особенности геометрии звеньев (рис. 4б). Затем, при q₁ = q₂ = ... = q_n = 0, находим направления осей СК звеньев (рис. 4в). В соответствии с рисунком определяем матрицы C_i и вектора l_i:

C₁ = C_z, C₂ = C₃ = C₄ = C_y, C₅ = C_x, C₆ = C_z;

l₁ =  ; l₂ = ; l₃ = ; l₄ =  ; l₅ =  ; l₆ = ;

(l_ix, l_iy, l_iz определяют геометрические размеры звеньев).

        Рисунок 4а.                          Рисунок 4б.                  Рисунок 4в.

Промышленный робот-манипулятор РМ-01 (или PUMA 560 в иностранной литературе) – антропоморфный манипулятор с 6 вращательными шарнирами.  На рис. 5 обозначены оси вращения его шарниров.

Рисунок 5.

Для описания кинематики робота РМ-01 свяжем с каждым из его звеньев соответствующие системы координат (рис. 6).

  Отметим, что СК S₀ – неподвижна, а ее начало совпадает с началом СК S₁, жестко связанной с первым звеном. На рис. 6 и 7 робот находится в положении q₁= 0^o, q₂= -90^o, q₃= 90^o, q₄= 0^o, q₅= 0^o, q₆= 0^o. В соответствии с рис. 6 и 7 определяем матрицы C_i и вектора l_i:

C₁ = C_z, C₂ = C₃ = C_y, C₄ = C_z , C₅ = C_y, C₆ = C_z;

l₁ =  ; l₂ = ; l₃ = ; l₄ =  ; l₅ =  ; l₆ = ;

(l_ix, l_iy, l_iz определяют геометрические размеры звеньев; l_6я – длина захватного устройства).

Рисунок 6.

 Далее рассмотрим применение описанного метода для получения основных соотношений кинематики манипуляторов.

Координаты точки G (схвата манипулятора) могут быть найдены из соотношения:

R_G = l_i⁰ =   C₁ l₁ + C₁ C₂ l₂  + … +  C₁ C₂ … C_n l_n = C_1i l_i         (2.1)

где принято обозначение C_1i = C₁ C₂ … C_i

l_i⁰  - координаты вектора l_i  в СК S₀  (см. рис. 8).

Рисунок 7.

Матрица C_G = C_1n определяет ориентацию СК схвата (совпадающую с СК S_n) относительно базовой СК S₀. Пара (R_G, C_G) задает решение прямой задачи кинематики манипулятора.

Линейная скорость схвата определяется выражением:

V_G = R_G¢ = C₁¢l₁ + (C₁¢C₂ + C₁C₂¢)l₂ +...+( C₁¢C_2...C_n+..+ C₁C_2…C_n¢)l_n      (2.2)

(здесь учтено, что l_i¢ = 0, т.к. эти вектора рассматриваются отно­сительно неподвижных систем координат).

Рисунок 8.

Производные от матриц поворотов можно вычислить с помощью вспомогательных матриц Q_i:

C_i¢ = Q_i C_i q_i¢

                   0 0  0      0 0 1       0 -1 0

где Q_i = { ( 0 0 -1 ), ( 0 0 0  ), ( 1  0 0  )} при C_i = {C_x, C_y, C_z}  соответственно

                   0 1  0      -1 0 0        0  0 0

Введем обозначение:

V_ps = C₁ C₂ ... C_s-1 Q_s C_s ... C_p

Тогда (2.2) примет вид:

V_G = V₁₁ q₁¢ l₁  +...+ (V_n1 q₁¢ +...+ V_nn q_n¢) l_n =( V_ps l_p) q_s¢        (2.3)

Угловая скорость схвата (конечного звена) по теореме сложения угловых скоростей равна геометрической сумме угловых скоростей звеньев:

W_G =  W_S =   e_{rot  s}⁰ q_s¢ =  C_1s e_{rot  s} q_s¢                               (2.4)

Теперь можно определить матрицу Якоби манипулятора, используемую при управлении для вычисления программных скоростей в шарнирах по желаемому движению в декартовом пространстве, а также для анализа вырожденных конфигураций манипулятора. По определению матрица Якоби манипулятора J(q) связывает вектора V_G и W_G с вектором обобщенных  координат  манипулятора q = (q₁, q₂, ..., q_n)^T   (символ “T” означает транспонирование):

(V_G, W_G ) = J(q) q¢                                                                                    (2.5)

Из (2.5), воспользовавшись выражениями (2.3), (2.4), можно опре­делить компоненты матрицы Якоби:

J_vs  =  V_ps l_p;  J_ws  = C_1s e_{rot  s}                                                        (2.6)

где   J_vs - компоненты   первых   трех   строк   матрицы    Якоби, J_ws  - компоненты последних трех строк.

В дальнейшем, при выводе уравнений динамики манипуляторов, нам понадобятся также выражения для радиус-вектора R_i и скорости V_i произвольной точки M_i i-ого звена манипулятора (см. рис. 9):

R_i = l_p⁰ + r_i⁰  = C_1p l_p + C_1i r_i

V_i = ( V_ps q_s¢) l_p  +  V_is q_s¢ r_i                   (r_i =О_iМ_i)                          (2.7)

Рисунок 9.

3. Формирование системы уравнений динамики

В этом параграфе будут получены уравнения динамики манипулятора с перпендикулярными или параллельными осями соседних шар­ниров. При их выводе будем использовать описание кинематики с помощью матриц поворотов и векторов переноса и уравнения Лагран­жа II рода:

 t_k                                                                           (3.1)

где L = (K - P) - функция Лагранжа, K и P - кинетическая и потенциальная энергия манипулятора; t_k - момент обобщенных сил в     k-ом шарнире, обусловленный работой привода и воздействием внешних нагрузок.

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

Для вычисления кинетической энергии воспользуемся формулой:

K = К_i,      К_i  =                                                                  (3.2)

где К_i - кинетическая энергия i-ого звена манипулятора; dК_i - ки­нетическая энергия элемента массы dm_i i-ого звена:

dК_i  = 1/2 V_i² dm_i  = 1/2 (V_i ,V_i ) dm_i                                                     (3.3)

Преобразуем выражение для V_i из соотношения (2.7):

V_i =  (V_ps q_s¢) l_p  +  V_is q_s¢ r_i  =  (V_ps q_s¢) r_p  

где r_p =      (i=1,..,n)

Обозначим w^p = V_ps q_s¢; тогда (3.3) можно представить в виде: 

dК_i =1/2 (w^p r_p , w^l r_l) dm_i = 1/2 Tr (w^p r_p r_l^Tw^{l T}) dm_i;

здесь Tr(A) - след матрицы А.

Тогда К_i  = 1/2 Tr (w^p X_pl w^{l T}),  где X_pl  =              (3.4)

причем X_pl  =                                                      (3.5)

 (здесь r_C,i  - центр масс i-ого звена в СК S_i;

 J_i      - матрица инерции i-ого звена манипулятора;

          см. приложение I)

С учетом (3.4) выражение для кинетической энергии манипулятора примет вид:

K = К_i  = 1/2 Tr (w^p X_pl w^{l T}). 

Теперь определим потенциальную энергию манипулятора как сумму потенциальных энергий его звеньев:

P =  P_i;  P_i = - m_i (g, r_C,i⁰  ) = - m_i (g,  C_1p l_p + C_1i r_C,i);

Найдем выражение для части уравнений (3.1), в которые вхо­дит кинетическая энергия манипулятора:

[Tr () + Tr () –

                                              Tr ()]

Заметим, что первое и третье слагаемые выражения в квадратных скобках равны между собой. Действительно,

 =                                             (3.6)

где обозначено .

                                                               (3.7)

Т.к. из определения величин V_pk   и V_pks  следует,  что  V_pks =  V_psk, то формулы (3.6) и (3.7) совпадают.

Итак, необходимо вычислить значение:

  Tr ()                                                                 (3.8)

Т.к. =V_pk, а == = + ,

то (3.8) можно представить в виде:

[ Tr () + Tr ()]

Подставив значения X_pl из формул (3.5) и сделав, где это необходимо, переход от функции Тr к скалярному произведению, получим:

+                     (3.9)

где r_p =    (i=1,..,n)                      

Потенциальная энергия манипулятора входит в уравнения динамики (3.1) в виде:

               (3.10)

Перегруппировав компоненты в формулах (3.9)-(3.10), можно полу­чить уравнения динамики в виде:

или в матричной форме:

                                                                      (3.11)

В уравнениях (3.11) D(q) - симметричная, положительно определенная матрица инерции манипулятора [60] с элементами:

d_ks  =                       (3.12)

- вектор кориолисовых и центробежных сил:

,

                   (3.13)

p - вектор гравитационных сил с компонентами:

                                                                    (3.14)

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

Приложение I. Определение матриц инерции звеньев манипулятора.

Матрица инерции i-ого звена манипулятора имеет вид:

J_i  = 

Компоненты этой матрицы можно выразить через компоненты (I_i)_ps тензора инерции i-ого звена [23]:

(J_i)_ps =  (I_i)_ps, при p¹s

(J_i)₁₁ = 1/2 [-(I_i)₁₁  + (I_i)₂₂  + (I_i)₃₃]

(J_i)₂₂ = 1/2 [ (I_i)₁₁   - (I_i)₂₂   + (I_i)₃₃]

(J_i)₃₃ = 1/2 [ (I_i)₁₁  + (I_i)₂₂  -  (I_i)₃₃]

4. Оценка вычислительной эффективности уравнений динамики

В [1, 5] и других работах, связанных с исследованием мето­дов формирования уравнений динамики манипуляторов, отмечается, что с вычислительной точки зрения уравнения в форме Лагранжа неэффективны по сравнению с другими способами описания динамики. При этом наиболее эффективные алгоритмы не представимы в замкну­том, матричном виде, что усложняет их использование в задачах управления роботами с учетом динамики.

Представленный в разделе 3 метод описания динамики, имеет мат­ричную структуру и обладает довольно высокой вычислительной эф­фективностью. Для ее оценки необходимо определить число арифметических операций, требуемых для реализации алгоритма, а также количество обращений к тригонометрическим функциям.

Распределение вычислительных затрат, необходимых при моделировании динамики манипулятора, представлено в таблице 2. Для  указано число операций, которые необходимо выполнить дополнительно после расчета d_ks. При оценивании объемов вычислений принимались во внимание известные свойства матриц D(q) и H_k = {}  [45]: d_ks  = d_sk; = 0 при k=s³t;  = ; = -     при k,s³t.

Отметим, что в суммарный объем вычислений входят также затраты на решение прямой задачи кинематики и расчет элементов матрицы Якоби манипулятора. Число обращений к тригонометрическим функциям - 2n.

Таблица 2.

Результаты, представленные в таблице 2, характеризуют объем вычислений для манипулятора общего вида. При рассмотрении манипулятора конкретного типа число операций сокращается, т.к. некоторые компоненты векторов l_i , r_{C, i} и матриц , V_ps и V_psk бу­дут иметь нулевые значения.

Для манипулятора БКМ вычислительные затраты на расчет полной модели динамики составляют 3479 опера­ций умножения и 2503 операции сложения. Программная реализация алгоритма расчета уравнений динамики позволяет учитывать особенности геометрии и распределения масс для конкретных манипулято­ров, что повышает вычислительную эффективность алгоритма.

5. Описание кинематики и динамики манипуляторов с поступательными шарнирами

Рассмотрим применение метода последовательного формирования систем координат звеньев с использованием матриц размера 3х3 для описания кинематики и динамики манипуляторов с поступательными шарнирами.

Выражения для пары (R_G ,С_G), которая определяет решение пря­мой кинематической задачи, будут иметь тот же вид, только теперь матрица преобразования систем координат соседних звеньев С_i бу­дет единичной матрицей, если i-ый шарнир - поступательный, а век­тора l_i не будут постоянными в СК S_i i-ого звена:

R_G  =   C_1i l_i;

C_G = C_1n;   C_i = {C_x, C_y, C_z} если i - вращательный шарнир

C_i º E                если i - поступательный шарнир

Введем параметр ind_s, характеризующий тип шарнира: ind_s =1 для вращательного шарнира, ind_s = 0 для поступательного. Тогда V_G можно представить в виде:

V_G = R_G¢ = C₁¢l₁ + C₁ l₁¢ ind₁ +...+( C₁¢_...C_n+..+ C_1…C_n¢)l_n + C₁C_2...C_nl_n¢ ind_n

Заметим, что l_p¢ = q _p¢ / q _p l _p , т.к. в силу выбора систем  координат звеньев перемещение  вдоль  поступательных  шарниров  всегда происходит вдоль одной  из  координатных  осей  соответствующего звена (l_p¢ = l_p =0 при q_p =0).

С учетом этого перепишем выражение для V_G:

V_G=(V₁₁ q₁¢ (1-ind₁)+C₁₁q₁¢ ind₁/q₁ )l₁+…+(V_n1 q₁¢ +..+V _nn q_n¢ +C_1n q_n ind_n/q_n)l_n

Введем в рассмотрение новые матрицы  , связанные с V_ps следующими соотношениями:

 =   

Тогда для V_G справедливы соотношения, аналогичные по виду (2.3):

 V_G =                                                                              (5.1)

Угловая скорость схвата при наличии поступательных шарниров имеет вид:

                                                      (5.2)

(5.1) и (5.2) определяют значения компонент матрицы Якоби:

          ;     

Прежде, чем получить соотношения для R_i и V_i - положения и скорости произвольной точки М_i i-ого звена манипулятора - отме­тим, что если i-ый шарнир поступательный, то введенный в разделе 2 вектор не будет постоянным в СК S_i , и, следовательно, выра­жения (2.7) в этом случае не могут быть использованы для формирования кинетической и потенциальной энергии манипулятора, т.к. матрицы J_i =  теряют физический смысл. Поэтому в ка­честве будем брать вектор  (см. рис. 10).

Рисунок 10.

Заметим, что если i-ый шарнир поступательный, то матрица инерции i-ого звена будет определяться в СК, полученной параллельным переносом S_i в точку О_i+1.

Радиус-вектор точки М_i представим в виде:

Тогда:

где r_p =     (i=1,..,n)

          Итак, для случая, когда манипулятор содержит поступательные шарниры, мы получили выражения для R_i и V_i в форме, аналогичной (2.7) из 2, когда предполагалось, что манипулятор содержит лишь вращательные шарниры. Поэтому все дальнейшие выкладки будут аналогичны рассмотренным выше при выводе уравнений динамики (отметим только, что в них необходимо использовать новые значения ). Следовательно, предлагаемый метод дает возмож­ность получить уравнения динамики и для манипуляторов с поступа­тельными шарнирами, причем вид этих уравнений совпадает с (3.11) - (3.14). Отметим, что объем вычислений, требуемых для реа­лизации алгоритма, в этом случае сократится, т.к. для поступа­тельных координат структура матриц С_ip и V_ps упрощается.

6. Применение символьных преобразований

Для повышения вычислительной эффективности уравнений динамики манипулятора были получены выражения для динамических коэф­фициентов в символьном виде с помощью языка аналитических вычис­лений REDUCE [40]. Разработан пакет программ, позволяющий полу­чать в символьном виде уравнения динамики произвольного манипу­лятора с заданным количеством степеней свободы.

Для получения символьных выражений используются соотношения (3.12)-(3.14). Входными параметрами программы на языке REDUCE являются кинематические, геометрические и масс-инерционные параметры манипулятора. Результатом работы программы являются урав­нения динамики манипулятора или выражения для коэффициентов мат­рицы инерции и моментов кориолисовых и центробежных сил инерции и силы тяжести в символьном виде. Имеется возможность генерации фрагментов программ расчета динамических параметров на языке ФОРТРАН.

Вычислительная эффективность получаемых соотношений анали­зировалась для моделей 2-х и 3-х степенных манипуляторов. Рассматривались манипуляторы, у которых вектора l_i  и r_Ci имеют одну ненулевую компоненту, а матрицы инерции звеньев диагональные с произвольными элементами.

Для 2-степенного манипулятора динамические коэффициенты имеют вид (приведен фрагмент выходного файла, сгенерированного программой символьной обработки):

Коэффициенты матрицы инерции манипулятора:

D(1,1)=cos(q2)*L1*L2*M2+I1Z+I2Z+L1**2*M2;

D(1,2)=(cos(q2)*L1*L2*M2+2*I2Z)/2;

D(2,1)=D(2,1); D(2,2)=I2Z;

Коэффициенты матриц кориолисовых и центробежных сил:

H(1,1,1)=0;

H(1,1,2)=-(sin(q2)*L1*L2*M2)/2;

H(1,2,1)=H(1,1,2);

H(1,2,2)=-(sin(q2)*L1*L2*M2)/2;

H(2,1,1)=-H(1,2,1); H(2,1,2)=0; H(2,2,1)=0; H(2,2,2)=0;

Компоненты вектора моментов силы тяжести:

P(1)=(G*(cos(q1)*cos(q2)*L2*M2+cos(q1)*L1*M1+

2*cos(q1)*L1*M2-sin(q1)*sin(q2)*L2*M2))/2; P(2)=(G*L2*M2*(cos(q1)*cos(q2)-sin(q1)*sin(q2)))/2;

Соответствующие выражения для 3-степенного манипулятора:

Коэффициенты матрицы инерции манипулятора:

D(1,1)=2*sin(q2)**2*sin(q3)**2*(-I3X+I3Z)+       sin(q2)**2*cos(q3)*L2*L3*M3+ sin(q2)**2*(I2X-I2Z+I3X-I3Z+L2**2*M3)+ 2*sin(q2)*sin(q3)*cos(q2)*cos(q3)*(I3X-I3Z)+ sin(q2)*sin(q3)*cos(q2)*L2*L3*M3+ sin(q3)**2*(I3X-I3Z)+I1Z+I2Z+I3Z;

D(1,2)=0; D(1,3)=0; D(2,1)=0;

D(2,2)=cos(q3)*L2*L3*M3+I2X+I3X+L2**2*M3;

D(2,3)=(cos(q3)*L2*L3*M3+2*I3X)/2;

D(3,1)=0; D(3,2)=D(2,3); D(3,3)=I3X;

Коэффициенты матриц кориолисовых и центробежных сил:

H(1,1,1)=0

H(1,1,2)=(4*sin(q2)**2*sin(q3)*cos(q3)*(-I3X+I3Z)-

2*sin(q2)**2*sin(q3)*L2*L3*M3+

4*sin(q2)*sin(q3)**2*cos(q2)*(-I3X+I3Z)+

2*sin(q2)*cos(q2)*cos(q3)*L2*L3*M3+

2*sin(q2)*cos(q2)*(I2X-I2Z+I3X-I3Z+L2**2*M3)+

2*sin(q3)*cos(q3)*(I3X-I3Z)+sin(q3)*L2*L3*M3)/2; H(1,1,3)=(4*sin(q2)**2*sin(q3)*cos(q3)*(-I3X+I3Z)-

sin(q2)**2*sin(q3)*L2*L3*M3+

4*sin(q2)*sin(q3)**2*cos(q2)*(-I3X+I3Z)+ sin(q2)*cos(q2)*cos(q3)*L2*L3*M3+ 2*sin(q2)*cos(q2)*(I3X-I3Z)+ 2*sin(q3)*cos(q3)*(I3X-I3Z))/2;

H(1,2,1)=H(1,1,2);

H(1,2,2)=0; H(1,2,3)=0; H(1,3,1)=H(1,1,3);

H(1,3,2)=0; H(1,3,3)=0; H(2,1,1)= -H(1,2,1);

H(2,1,2)=0; H(2,1,3)=0; H(2,2,1)=0; H(2,2,2)=0;

H(2,2,3)= -(sin(q3)*L2*L3*M3)/2;

H(2,3,1)=0; H(2,3,2)=H(2,2,3); H(2,3,3)= -(sin(q3)*L2*L3*M3)/2;

H(3,1,1)= -H(1,3,1); H(3,1,2)=0; H(3,1,3)=0;

H(3,2,1)=0; H(3,2,2)= -H(2,3,2);

H(3,2,3)=0; H(3,3,1)=0; H(3,3,2)=0; H(3,3,3)=0;

Компоненты вектора моментов силы тяжести: P(1)=(sin(q2)*cos(q1)*cos(q3)*G*L3*M3+

sin(q2)*cos(q1)*G*L2*(M2+2*M3)+

sin(q3)*cos(q1)*cos(q2)*G*L3*M3)/2;

P(2)=(-sin(q1)*sin(q2)*sin(q3)*G*L3*M3+

sin(q1)*cos(q2)*cos(q3)*G*L3*M3+

sin(q1)*cos(q2)*G*L2*(M2+2*M3))/2;

P(3)=(-sin(q1)*sin(q2)*sin(q3)*G*L3*M3+

sin(q1)*cos(q2)*cos(q3)*G*L3*M3)/2;

Для 2-степенного манипулятора рассматриваемого класса объем вычислительных затрат составляет 135 и 76 сложений. Для приведенных символьных выражений: 34 и 9, т.о. вычислительная эффективность повысилась в 5 раз. Для 3-степенного манипулятора имеем соответственно 412 и 259 для численной модели и 135 и 47 для символьной; повышение эффективности в 4 раза.

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

7. Модели приводов и механических передач

Для замыкания системы уравнений (3.11) необходимо получить выражения для обобщенных моментов в шарнирах. Они определяются типами и параметрами двигателей, механических передач (редукто­ров), а также особенностями системы управления робота. Здесь будут рассмотрены электромеханические привода с двигателями пос­тоянного тока, обратимыми зубчатыми редукторами и замкнутыми по скорости следящими системами. Будут представлены различные по сложности модели приводов, в которых учитываются упругость в шарнирах и нелинейные элементы: люфт, сухое трение, муфты пре­дельного момента, тормоза.

Рассмотрим i-ый шарнир манипулятора. Баланс моментов для него имеет вид:

                                                (7.1)

В этих уравнениях  и - момент инерции якоря двигателя и развиваемый двигателем электромагнитный момент, приведенные к выходу редуктора:  = , , где - передаточное число редуктора, равное отношению угловых скоростей звена и двигателя,  - коэффициент пропорциональности момента,  - ток в дви­гателе. Коэффициенты матрицы инерции манипулятора , компоненты вектора кориолисовых и центробежных сил , и гравитационных сил рассчитываются по формулам (3.12)-(3.14) соответственно.

Уравнения для двигателя постоянного тока имеют вид:

где  - соответственно напряжение, сопротивление и индуктивность в обмотках двигателя;  - коэффициент противо­-э.д.с., . Кроме того, для скоростной следящей системы:

,

где  - программная скорость i-ого звена.

Итак, имеем следующую систему уравнений, описывающих движе­ние i-ого звена:

Тогда для манипулятора справедливы векторные уравнения:

                                                                    (7.4)

Диагональные элементы матрицы  равны ; , ,  и  - диагональные матрицы (n x n) с элементами , , и  соответственно.

          Для позиционно-скоростной следящей системы имеем:

,

где  - программное положение i-ого звена.

Тогда справедливы векторные уравнения:

                                                 (7.4¢)

Здесь  - диагональная матрица (n x n) с элементами .

Системы (7.4) и (7.4¢) содержат уравнения третьего порядка, однако во многих случаях (см., например, [5]) достаточно рассматривать модели второго порядка для приводов, т.е. пренебрегать членами , т.к. время переходных процессов, связанных с перерегулиро­ванием тока, мало по сравнению с длительностью механических пе­реходных процессов вследствие значительной инерционности манипу­ляционной системы.

Если механические передачи манипулятора обладают люфтом, то для однозначного описания динамики манипулятора в этом случае необходимо ввести дополнительные переменные  и , задающие положение и угловую скорость валов двигателей, которые вместе с  и W  определяют в любой момент состояние манипулятора ( и  приведены к выходу редукторов). Обозначим через  половину величины люфта на выходе передачи. Тогда для всех i=1,..,n в случае, если люфт не выбран, т.е. при , движение в i-ом шарнире задается уравнениями:

                                                        (7.5)

В противном случае, движение будет определяться уравнениями (7.1). Кроме того, для шарниров, в которых выбран люфт, в каждый момент времени необходимо следить за условием сохранения связи ( -момент сил реакции на выходе передаточного механизма), и, в случае его нарушения, переходить к системе (7.5).

Рассмотрим теперь движение манипулятора при наличии моментов сопротивления , создаваемых силами сухого трения на выход­ных валах механических передач (учет сил вязкого трения не пред­ставляет труда, т.к. они входят в уравнения в виде добавочного коэффициента к пропорциональным обобщенным скоростям членам). Система уравнений динамики манипулятора примет вид:

                   (7.6)

При  момент  вычисляется по формуле , где - величина предельного момента сил сухого трения в i-ой передаче.

Если в процессе движения станет , то для вычисления  нужно сравнить  (модуль момента всех сил, действующих в i-ом шарнире), и . Если , то ; в  противном случае по i-ой координате будет выполняться условие  до мо­мента, пока не станет .

При наличии упругости в шарнирах уравнения (7.1) примут вид:

где - суммарный коэффициент упругости  i-ого  шарнира.  Если  в этом случае учесть также люфт в редукторе, то получим систему:

                                             (7.7)

                                                      (7.8)

Заметим, что, в отличие от рассмотренной выше модели привода с люфтом без учета упругости, здесь не требуется следить за вели­чиной реакции в передаче для определения момента разрыва связи, т.к. его значение может быть найдено из анализа двух независимых переменных  и : при связь окажется нарушенной.

        Рассмотрим теперь модель динамики манипулятора, который содержит на выходе редукторов муфты предельного момента, имеющие линейную характеристику с насыщением [47] (рис. 11). Здесь - значение момента  в передаче, когда  в  муфте   начнется проскальзывание: если , то  и движение манипулятора  описывается  системой (7.7);  в  противном  случае .

Рисунок 11.

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

         для

Движение i-ого звена манипулятора можно определить из уравнения:

где  рассчитывается  по  формулам  (7.8)  с  учетом  условия .

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

Выводы.

Предложено математическое описание алгоритма последовательного формирования локальных систем координат звеньев для манипу­ляторов с вращательными и поступательными шарнирами, соседние оси которых параллельны или перпендикулярны. Разработан метод формирования уравнений динамики для такого класса манипуляторов с помощью уравнений Лагранжа II рода.

Уравнения имеют матричный вид, позволяют решать прямую и обратную задачи динамики, удобны для реализации на ЭВМ. Использование матриц размера 3х3 обеспечило высокую вычислительную эффективность уравнений; разработанная программная реализация алгоритма расчета коэффициентов уравнений динамики позволяет проводить оптимизацию вычислительных затрат для конкретных типов манипуляторов.

Структура полученных уравнений допускает применение к ним символьных преобразований. Разработан пакет программ на языке REDUCE для вывода уравнений динамики в символьном виде. Их вы­числительная эффективность в 4-5 раз превышает эффективность численных моделей для 2-х и 3-х звенных манипуляторов.

Отметим, что, несмотря на постоянный рост мощности компьютеров, требование высокой вычислительной эффективности уравнений динамики остается критичным. Это объясняется тем, что, во-первых, в системах управления роботов используются как правило относительно медленные процессоры, и для решения уравнений динамики в реальном времени необходимы эффективные алгоритмы расчета. А, во-вторых, сложность механических структур современных роботов (параллельных, с избыточными степенями подвижности, так называемых роботов-гуманоидов) требует эффективных методов расчета динамики для задач их моделирования и управления.

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

Литература

1.     Фу К., Гонсалес Р., Ли К. Робототехника.- М.: Мир, 1989.

2.     Kane T., Dynamics, New York, Holt,  Rihehart  and  Wiston, 1968.

3.     Виттенбург Й. Динамика систем твердых тел.- М.: Мир,1980.

4.     Denavit J, Hartenberg R.S. A kinematic notation for lower­-pair mechanisms based on matrices.,  J. Appl. Mech., 77, 1955, c.215-221.

5.     Вукобратович М., Стокич Д., Кирчански Н. Неадаптивное и адаптивное управление манипуляционными роботами. - М.: Мир, 1989.

6.     Попов Е.П. Управление роботами-манипуляторами. Изв.АН СССР, Техн. киберн., 1974, N 6, с.51-56.

7.     Vukobratovic M., Stepanenko Y. Mathematical model of gene­ral anthropomorphic systems. Math Biosciences, Vol.17, 1973, c.191-242.

8.     Накано Э. Введение в робототехнику, М.: Мир, 1988.

9.     Hollerbach J. A recursive Lagrangian formulation of manipu­lator dynamics and comparative study of dynamic complication complexity. IEEE Trans. on SMC, SMC-10, No 11, 1980, c.730-736.

10. Kahn M.E., Roth B. The near-minimum-time control of open -loop articulated kinematic chains, ASME J. of  Dynam Syst, Measur.and Countr., vol. 93, 1971, c.164-172.

11. Uicer J.J. Dynamic force analysis of spatial linkages, ASME J. of appl. mech., June, 1967, c.418-424.

12. Lee C.S.G., Lee B.H., Nigam R. Development of generalized d'Alambert Equation of motion for mechanical manipulators, Proc 2nd conf. Decision and Control, San Antonio, 1983, c. 1205-1210

13. Thomas M, Tesar D. Dynamic modeling of serial manipulator arms. Trans. of ASME, vol. 104, Sept, 1982,c.218-228.

14. Mahil S. On the application of Lagrange's method to the description of dynamic systems. IEEE Trans. on SMC, vol SMC-12, N 6, 1982.

15. Wang L.T., Ravani B. Recursive computations of kinematic and dynamic equations for mechanical manipulators. IEEE J. of Rob. and Autom., vol. RA-1, N 3, Sept. 1985, c.124-131.

16. Balafoutis C, Patel R., Misra P. Efficient modeling and computation of manipulator dynamics using orthogonal cartesian tensors. IEEE J. of Rob. and Autom., 4, N 6, c.665-676.

17. Castelain J.M, Bernier D. A new program based on the hiper­complex theory for automatic generation of the direct differen­tial model of robot manipulators . Mech. mach. theory, vol. 25, N 1, 1990, c.69-83.

18. Mladenova C. Mathematical modeling and control of manipu­lator systems. Int. J. Robotics and computer-integrated manufac­turing, vol. 8, N 4, 1991, c 233-242.

19. F.C. Park, J. Choi, and S.R. Ploen, ”A Li Group Formulation of Robot Dynamics,”The Int. J. of Robotics Research, Vol.14, No.6, Dec.1995.

20. Paul R. Manipulator cartesian path control. IEEE Trans.  on SMC-9, Febr, 1979, c.702-711.

21. Vukobratovic M, Potkonjak V. Contribution to automatic forming of active chain models via Lagrangian form. J of Appl. Mech., N 1, 1979.

22. Renaud N. An efficient iterative analytical procedure for obtaining a robot manipulator dynamic model. Proc. of FirstInt. Symp. of Rob. Research, Bretton Woods, New Hampshire, USA,1983.

23. Li C.G. A new method for dynamic analysis of robot manipu­lators . IEEE Trans. on Syst., Man and Cybern., 1988, 18, N 1, c.105-114.

24. Walker M.W., Orin D.E. Efficient dynamic computer simula­tion of robotic mechanisms. ASME J. of Dyn. Syst.,Meas. and Contr., vol. 104, Sept. 1982, c.205-211.

25. Armstrong W.W. Recursive solution to the equations of motion of an n-link manipulator. Proc of the 5th World Congress on Theory of Mach. and Mech, Montreal, 1979, c. 1343-1346.

26. Малышев А.Б., Чуменко В.Н. Универсальные программы модели­рования динамики манипуляционного робота. "Роботы и РТС", Иркутск, 1983, 117-126.

27. Попов Е.П., Верещагин А.Ф., Зенкевич С.Л. Манипуляционные роботы: динамика и алгоритмы.- М.: Наука, 1980.

28. Huston R.L. The use of Kane's metod in the modeling and simulation of robotic systems. Proc. IMACS Symp. Syst. Modeling and Simul., Cetraro, 18-21 sept, 1988.

29. Ma X., Xu X. A futher study of Kane's equations. Proc IEEE Int Conf Syst, Man and Cybern, Beijing, Aug. 8-12, 1988, c.107-112.

30. Коноплев В.А. Агрегативные модели механики систем твердых тел со структурой дерева. Изв. АН СССР, МТТ, N 6, 1989, с 46-54.

31. Погорелов Д.Ю., "Алгоритмы синтеза и численного интегрирования уравнений движения систем тел с большим числом степеней свободы", VIII Всероссийский съезд по теоретической и прикладной механике, Пермь, 2001, с. 490.

32. И.Ю.Балабан, Г.К.Боровин, В.В.Сазонов, “Язык программирования правых частей уравнений движения сложных механических систем”, Препринт ИПМ им. М.В. Келдыша РАН, N 62, 1998, 22 с.

33. Dapper, R. Maafl, V. Zahn, R. Eckmiller, Neural Force Control (NFC) Applied to Industrial Manipulators in Interaction with Moving Rigid Objects, Proceedings of the 1998 IEEE International Conference on Robotics & Automation, Leuven, Belgium,  May 1998.

34. S. Jung, S. B. Yim, T. C. Hsia, Experimental Studies of Neural Network Impedance Force Control for Robot Manipulators, Proceedings of the 2001 IEEE International Conference on Robotics & Automation, Seoul, Korea, May 2001.

35. G. Rodriguez, A. Jain and K. Kreutz-Delgado, "A Spatial Operator Algebra for Manipulator Modelling and Control," Int. J. Robotics Research, vol. 10, no. 4, pp. 371-381, 1991.

36. A. Jain, G. Rodriguez, Computational Robot Dynamics Using Spatial Operators, Proceedings of the 2000 IEEE International Conference on Robotics & Automation, San Francisco, CA, April 2000.

37. M. Emami, A. Goldenberg, I. Turksen, Fuzzy-Logic Dynamics Modeling of Robot Manipulators, Proceedings of the 1998 IEEE International Conference on Robotics & Automation, Leuven, Belgium,  May 1998.

38. R. Featherstone, D. Orin, Robot Dynamics: Equations and Algorithms, Proceedings of the 2000 IEEE International Conference on Robotics & Automation, San Francisco, CA, April 2000.

39. Cheng P., Weng C., Chen C. Symbolic derivation of dynamic equation of motion for robot manipulator using program symbolic method. IEEE J. Rob. and Autom, 4, N 6, 1988, c. 599-609.

40. Ju M.S., Mansor J.M. Comparision of methods for developing the dynamics of rigid body systems. Int. J. Rob. Res., N6, 1989, c.19-27.

41. Белоусов И.Р., “Применение метода символьных преобразований для формирования алгоритмов параллельных вычислений в задачах кинематики и динамики роботов”, Отчет ИПМ им. М.В. Келдыша РАН № 5-19-93, 1993, 25 с.

42. Lathrop L.H. Parallelism in manipulator dynamics.  Int.  J. Rob. Res., vol.4, No 2, 1985, c.80-102.

43. R. Featherstone, "A Divide-and-Conquer Articulated-Body Algorithm for Parallel O(log(n)) Calculation of Rigid-Body Dynamics. Part 1: Basic Algorithm," Int. Y. Robotics Research, vol. 18, no. 9, pp. 867-875, 1999.

44. A. Fijany, I. Sharf and G. M. T. D'Eleuterio, "Parallel O(logN) Algorithms for Computation of Manipulator Forward Dynamics," IEEE Trans. Robotics & Automation, vol. 11, no. 3, pp. 389-400, June 1995.

45. Vukobratovic M, Kircanski N, Real-time dynamics of manipu­lation robots, Springer-Verlag, 1985.

46. Han J.-Y. Fault-tolerant computing for robot kinematics using linear arithmetic code. IEEE Int. Conf. Robotics and Auto­mation, Cincinnati, May May 13-18, 1990, vol. 1, c.285-290.

47. Справочник по промышленной робототехнике.- М.: Машиностроение, 1990.