Привет, Хабр! Одна из задач при управлении роботами-манипуляторами – расчет обратной кинематики. Данный вид кинематики позволяет вычислить углы наклона суставов робота (joints) таким образом, чтобы захват (grip) робота пришел в заданные трехмерные координаты с правильным углом наклона. Для многих роботов уже есть алгоритмы и формулы вычисления обратной кинематики, мы (команда Zebrains) столкнулись с отсутствием готового решения для робота xArm 2.0.

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

Внешний вид робота xArm 2.0
Внешний вид робота xArm 2.0

Начнем с особенностей робота:

Согласно описанию, он имеет 6 степеней свободы, одна из которых относится к подвижной клешне. В нашем проекте было необходимо подбирать объекты различных форм, которые двигались по конвейеру и для сбора объектов клешня была заменена на присоску с насосом для вакуумного захвата. Таким образом, степеней свободны стало 5. Последний сустав робота отвечает за поворот присоски влево и вправо, что не требуется при вертикальном подборе объектов. В итоге, остается 4 степени свободы, управляя которыми необходимо перемещать захват в нужное положение.

Внешний вид робота xArm 2.0
Внешний вид робота xArm 2.0

Сначала мы пытались использовать различные готовые алгоритмы вычисления обратной кинематики, представленные в MoveIt2. Но местные алгоритмы предназначены для 6-ти осевых роботов (без учета поворота захвата), да при том, с несколько иным строением суставов (см. робот Panda на картинке ниже).

Внешний вид робота Panda
Внешний вид робота Panda

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

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

Демонстрация движения робота в rviz
Демонстрация движения робота в rviz

Решение задачи кинематики робота-манипулятора (кинематика №1)

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

Геометрическая схема расчетов 1-ой кинематики
Геометрическая схема расчетов 1-ой кинематики

При вычислении задачи обратной кинематики использовались следующие условия:

За поворот в сторону цели отвечает нулевой джоинт, поэтому для него угол поворота считаем отдельно, а дальше работаем в двухмерном пространстве;

Принимаем, что положение второго джоинта и цели, куда робот должен дотянуться лежат на одной прямой (P_2, P_5), а звено (P_3, P_4) параллельно этой прямой. Угол поворота прямой P_2,P_5) — ∠γ считается отдельно и влияет на угол поворота 2 джоинта.

Угол поворота нулевого джоинта вычисляется по формуле:

\angle 0 = \arctan\!\left(\frac{goal_y - robot_y}{goal_x - robot_x}\right)

где  \mathrm{goal}_x и \mathrm{goal}_y  — координаты точки назначения по x и y

\mathrm{robot}_x и\mathrm{robot}_y — координаты робота назначения по x и y

Поскольку робот может иметь свое отклонение по данному джоинту, то дополнительно добавим коэффициент смещения угла \angle 0 \mathrm{angle}_{0_{\mathrm{add}}} , а также коэффициент направления —  \mathrm{angle}_{0_{\mathrm{direction}}} \; \left( \mathrm{angle}_{0_{\mathrm{direction}}} \in \{1; -1\} \right).

Итоговая формула будет иметь следующий вид:

 \angle_0 = \mathrm{angle}_{0_{\mathrm{direction}}} \cdot \arctan \left( \frac{ \mathrm{goal}_y - \mathrm{robot}_y }{ \mathrm{goal}_x - \mathrm{robot}_x } \right) + \mathrm{angle}_{0_{\mathrm{add}}}.

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

 P_5(x) = \sqrt{ \left( \mathrm{robot}_x - \mathrm{goal}_x \right)^2 + \left( \mathrm{robot}_y - \mathrm{goal}_y \right)^2 }  P_5(y) = \mathrm{goal}_z - d_0 - \mathrm{robot}_z

Где d_0— длина нулевого джоинта,

goal_z — координата точки назначения по оси z,

robot_z — координата робота по оси z.

Расстояние dist рассчитывается по формуле:

 \text{dist} = \sqrt{ P_5(x)^2 + P_5(y)^2 }

Примем новую систему координат, в которой dist лежит на оси x, а перпендикуляры к этой прямой на оси y.

Тогда, точки P_3 и P_4 лежат на одной прямой и имеют высоту h=y. Эта высота присутствует в треугольниках △P_2P_3A  и △P_4P_5B. Через эти же треугольники можно найти углы ∠α_2(поскольку он равен углу ∠P_3P_2A, как углы, образованные общей прямой (P_2,P_3) и параллельными отрезками (P_2,A) и (P_3,P_4) и ∠α_3 (равен углу ∠P_4P_5B как накрест лежащие, образованные параллельными отрезками (B,P_5)и (P_4,N_3) и пересекающим их отрезком (P_4,P_5).

Для нахождения углов ∠α_1 и ∠α_3  необходимо знать значения катетов P_2,A =b_1иB,P_5=b_2,гипотенузы (P_2,P_3)=a_1 и (P_4,P_5)=a_3 уже известны — это длины звеньев робота (между 2 и 3 джоинтами, между 4 и 5 джоинтами). Также, через уравнения окружностей (приняв P2за начало координат):

Уравнение окружности с центром в начале координат имеет вид:

a_1^2 = x_1^2 + y_1^2a_3^2 = (x_2 - \text{dist})^2 + y_2^2

Где  x_1, y_1 — координаты точки P_3

x_2, y_2— координаты точки P_4  а так как выше мы установили, что они лежат на одной прямой, то y_1=y_2,x_1=b_1,x_2=b_1+a_2

Следовательно:

a_1^2 = b_1^2 + y^2a_3^2 = (b_1 + a_2 - \text{dist})^2 + y^2

Поскольку y_2есть в обоих уравнениях, приравняем уравнения, выразив их через y_2

a_1^2-b_1^2=a_3^2-(b_1+a_2-dist )^2a_1^2-b_1^2=a_3^2-(b_1^2+2⋅b_1⋅(a_2-dist )+(a_2-dist)^2 )a_1^2-b_1^2=a_3^2-b_1^2-2⋅b_1⋅(a_2-dist)-(a_2-dist )^22⋅b_1⋅(a_2-dist )=a_3^2-a_1^2-(a_2-dist )^2b_1=(a_3^2-a_1^2-(a_2-dist )^2)/(2⋅(a_2-dist ) ).

Тогда:

b_2=dist-b_1-a_2.

Вычисления углов:

\angle \alpha_2 = 90^\circ - \arcsin\left(\frac{b_1}{a_1}\right)\angle \alpha_3 = \arccos\left(\frac{b_2}{a_3}\right)\angle \ y = \arcsin\left(\frac{y}{dist}\right)\angle \ a_1 = \arcsin\left(\frac{b_1}{a_1}\right)-∠y,

Где y - высота до цели относительно второго джоинта в мировых координатах.

Финальные формулы для вычисления углов поворота джоинтов с учетом коэффициента смещения имеют следующий вид:

\angle1 = \text{angle}_{1\,\text{direction}} \cdot \left( \arcsin\left(\frac{b_1}{a_1}\right) - \angle \ y \right)\angle2 = \text{angle}_{2\,\text{direction}} \cdot \left( 90^\circ - \arcsin\left(\frac{b_1}{a_1}\right) \right) \angle3 = \text{angle}_{3\,\text{direction}} \cdot \left( \arccos\left(\frac{b_2}{a_3}\right) \right)

Фрагмент кода, отвечающий за просчет кинематики №1:

bool calcIKSimple(float goalX, float goalY, float goalZ, const IkParameters &ikParams,
                  std::vector<double> &jointGroupPositions)
{
    double localX = sqrt(pow((ikParams.robotX - goalX), 2) + pow((ikParams.robotY - goalY), 2));
    double localY = goalZ - ikParams.d0 - ikParams.robotZ;
    double distance = sqrt(pow((localX), 2) + pow((localY), 2));

    jointGroupPositions[ikParams.joints[0]] =
            ikParams.angleZeroDirection * atan2(goalY - ikParams.robotY, goalX - ikParams.robotX)
            + ikParams.angleZeroAdd;
    jointGroupPositions[ikParams.joints[1]] = 0.0;
    jointGroupPositions[ikParams.joints[2]] = 0.0;
    jointGroupPositions[ikParams.joints[3]] = 0.0;
    jointGroupPositions[ikParams.joints[4]] = 0.0;

    if (distance > ikParams.d3 && distance < ikParams.d1 + ikParams.d2 + ikParams.d3) {
        double b1 = (pow(ikParams.d1, 2) - pow(ikParams.d3, 2) + pow((distance - ikParams.d2), 2))
                / (2 * (distance - ikParams.d2));
        double b3 = distance - b1 - ikParams.d2;

        double gamma = std::asin(localY / distance);
        jointGroupPositions[ikParams.joints[1]] = std::asin(b1 / ikParams.d1);
        jointGroupPositions[ikParams.joints[2]] =
                ikParams.angleTwoDirection * (M_PI_2 - jointGroupPositions[ikParams.joints[1]]);
        jointGroupPositions[ikParams.joints[1]] =
                ikParams.angleOneDirection * (jointGroupPositions[ikParams.joints[1]] - gamma);

        jointGroupPositions[ikParams.joints[3]] =
                ikParams.angleThreeDirection * std::acos(b3 / ikParams.d3);
        return true;
    } else {
        return false;
    }
}

Решение задачи кинематики робота-манипулятора (кинематика №2)

Назначение второго вида кинематики — расчет углов поворота джоинтов таким образом, чтобы робот пришел в заданную точку с соблюдением заданного угла последнего джоинта относительно горизонтальной линии.

Заданный угол — ∠β.

Геометрическая схема расчетов 2-ой кинематики представлена на рисунке 6

Рисунок 6 –  Геометрическая схема расчетов 2-ой кинематики
Рисунок 6 –  Геометрическая схема расчетов 2-ой кинематики

Найдем координаты точки P_4. Длина отрезка (P_3, P_4) равна длине джоинта 3 — a_3. Точка P_5 — цель, куда должен приехать робот. Ее координаты мы знаем. Рассмотрим треугольник △P_4P_5N_1. В нем угол △P_4P_5N_1 известен и равен ∠β. Таким образом, координаты точки P_4P_5N_1 можно получить по следующим формулам:

P_4 (x)=P_5 (x)-cos(∠β)⋅a_3;P_4 (y)=sin(∠β)⋅a_3+P_5 (y).

Рассмотрим треугольник △P_2P_3P_4. В нем нам известны длины всех трех сторон:

(P_2,P_3 )=a_1;(P_3,P_4 )=a_2;(P_2,P_4 )=√(P_4 (x)^2+P_4 (y)^2 )=d_1.

Чтобы найти угол отклонения первого джоинта ∠a_1 необходимо найти углы ∠y_3и ∠y_2, так как:

90^∘=∠α_1+∠γ_3+∠γ_2;∠α_1=90^∘-∠γ_3-∠γ_2.

По теореме косинусов [ [8], [9], [10] ] найдем угол ∠y_3

\angle \gamma_3 = \arccos\left(\frac{d_1^2 + a_1^2 - a_2^2}{2 \cdot d_1 \cdot a_1}\right)

:

Рассмотрим треугольник △P_2P_4N_2. В нем:

(P_2,P_4 )=d_1;(P_4,N_2 )=P_4 (y).

Таким образом, угол ∠y_2 можно рассчитать по формуле:

\angle \gamma_2 = \arcsin\left(\frac{P_4(y)}{d_1}\right)

Наконец, получаем формулу вычисления угла первого джоинта:

\angle \alpha_1 = 90^\circ - \arccos\left(\frac{d_1^2 + a_1^2 - a_2^2}{2 \cdot d_1 \cdot a_1}\right) - \arcsin\left(\frac{P_4(y)}{d_1}\right)

Финальная формула для вычисления угла поворота первого джоинта с учетом коэффициента смещения имеет следующий вид:

\angle1 = \text{angle}_{1\,\text{direction}} \cdot \left( 90^\circ - \arccos\left(\frac{d_1^2 + a_1^2 - a_2^2}{2 \cdot d_1 \cdot a_1}\right) - \arcsin\left(\frac{P_4(y)}{d_1}\right) \right)

Для вычисления угла отклонения второго джоинта ∠a_2потребуется вычислить угол ∠y_4, который является смежным углу ∠a_2, так как угол ∠P_4P_3N_5 — развернутый и равен 180°, то угол ∠a_2 можно найти по формуле:

180^∘=∠γ_4+∠α_2;∠α_2=180^∘-∠γ_4.

Повторно рассмотрим треугольник ∠(P_2,P_4,N_2 ) . По теореме косинусов получаем:

\angle \gamma_4 = \arccos\left(\frac{a_1^2 + a_2^2 - d_1^2}{2 \cdot a_1 \cdot a_2}\right)

Угол ∠a_2 вычисляется по формуле:

\angle \alpha_2 = 180^\circ - \arccos\left(\frac{a_1^2 + a_2^2 - d_1^2}{2 \cdot a_1 \cdot a_2}\right)

Финальная формула для вычисления угла поворота второго джоинта с учетом коэффициента смещения имеет следующий вид:

\angle2 = \text{angle}_{2\,\text{direction}} \cdot \left( 180^\circ - \arccos\left(\frac{a_1^2 + a_2^2 - d_1^2}{2 \cdot a_1 \cdot a_2}\right) \right)

Для нахождения третьего угла необходимо добавить несколько проекций:

Из точки P_3построим вертикаль. Поскольку отрезок (P_3,N_5) является продолжением отрезка (P_2,P_3), то угол ∠(N_3,P_3,N_5) равен углу ∠a_1

Также построим вертикаль из точки N_4 и продлим отрезок (P_4,N_4).

Таким образом угол ∠y_5будет равен сумме углов ∠a_1 и ∠a_2.

Рассмотрим треугольник P_4N_4P_5. Углы ∠y_5и ∠P_4N_4P_5 — вертикальные, а значит равны. Угол ∠y_1 вычисляется по формуле:

∠γ_1=90^∘-∠β.

Угол ∠a_3  вычислим зная, что сумма углов треугольника равна 180°, а два других угла треугольника нам известны. Получим:

∠α_3=180^∘-(∠α_1+∠α_2 )-(90^∘-∠β);∠α_3=90^∘-∠α_1-∠α_2+∠β.

Тут важно обратить внимание, что используются углы ∠a_1 и ∠a_2, т.е. углы поворота первого и второго джоинтов без учета направления. Финальная формула для вычисления угла поворота третьего джоинта с учетом коэффициента смещения имеет следующий вид:

\angle3= \text{angle}_{3\,\text{direction}} \cdot \left( 90^\circ - \angle \alpha_1 - \angle \alpha_2 + \angle \beta \right)

Фрагмент кода, отвечающий за просчет кинематики №2

bool calcIKFinalAngle(float goalX, float goalY, float goalZ, float angle,
    const IkParameters &ikParams, std::vector<double> &jointGroupPositions)
{
    double localX = sqrt(pow((ikParams.robotX - goalX), 2) + pow((ikParams.robotY - goalY), 2));
    double localY = goalZ - ikParams.d0 - ikParams.robotZ;

    jointGroupPositions[ikParams.joints[0]] =
        ikParams.angleZeroDirection * atan2(goalY - ikParams.robotY, goalX - ikParams.robotX)
        + ikParams.angleZeroAdd;
    jointGroupPositions[ikParams.joints[1]] = 0.0;
    jointGroupPositions[ikParams.joints[2]] = 0.0;
    jointGroupPositions[ikParams.joints[3]] = 0.0;
    jointGroupPositions[ikParams.joints[4]] = 0.0;

    double yP4 = std::sin(angle) * ikParams.d3 + localY;
    double xP4 = localX - std::cos(angle) * ikParams.d3;
    double distanceSmall = sqrt(pow((xP4), 2) + pow((yP4), 2));

    if (distanceSmall < ikParams.d1 + ikParams.d2) {

        jointGroupPositions[ikParams.joints[1]] = M_PI_2
            - std::acos((pow(distanceSmall, 2) + pow(ikParams.d1, 2) - pow(ikParams.d2, 2))
                        / (2 * distanceSmall * ikParams.d1))
            - std::asin(yP4 / distanceSmall);
        jointGroupPositions[ikParams.joints[2]] = M_PI
            - std::acos((pow(ikParams.d1, 2) + pow(ikParams.d2, 2) - pow(distanceSmall, 2))
                        / (2 * ikParams.d1 * ikParams.d2));
        jointGroupPositions[ikParams.joints[3]] = ikParams.angleThreeDirection
            * (M_PI - (M_PI_2 - angle)
               - (jointGroupPositions[ikParams.joints[1]]
                  + jointGroupPositions[ikParams.joints[2]]));

        jointGroupPositions[ikParams.joints[1]] =
            ikParams.angleOneDirection * jointGroupPositions[ikParams.joints[1]];
        jointGroupPositions[ikParams.joints[2]] =
            ikParams.angleTwoDirection * jointGroupPositions[ikParams.joints[2]];
        return true;
    } else {
        return false;
    }
}

Расчет рабочей области робота

Поскольку соблюдение заданного угла при перемещении в точку, накладывает на робота ограничение, следовательно, его рабочая зона сокращается с суммы длин его 2, 3 и 4 звеньев до меньшего значения, которое необходимо посчитать. Для этого, необходимо получить расстояние D =(P_2,N_2 )+(N_1,P_5 ). Точка P_5будет иметь наибольшее значение x при том же y и угле ∠β, если угол ∠α_2 будет равен 0, т.е. звенья a_1 и a_2 будут лежать на одной прямой, а (P_2,P_4 )=a_1+a_2. Тогда, в треугольнике △P_2 P_4 N_2 неизвестен только нужный нам катет (P_2,N_2 ), который можно найти следующим образом:

(P_2,N_2 )=√((P_2,P_4 )^2-P_4 (y)^2 ).

Из расчетов ранее, известно, что:

P_4 (y)=sin(∠β)⋅a_3+P_5 (y).

Рассмотрим треугольник P_4 P_5 N_1:

(P_4,N_1 )=P_4 (y)-P_5 (y);(P_4,P_5 )=a_3.

Тогда неизвестный катет находим по формуле:

(N_1,P_5 )=√((P_4,P_5 )^2-(P_4,N_1 )^2 );(N_1,P_5 )=√(a_3^2-(P_4 (y)-P_5 (y))^2 ).

Итоговая формула будет иметь следующий вид:

D=√((P_2,P_4 )^2-P_4 (y)^2 )+√(a_3^2-(P_4 (y)-P_5 (y))^2 ).

Результат разработки

Полученная кинематика была сначала протестирована в 3D игровом движке Unreal Engine 4, поскольку в нем при помощи системы Blueprints (визуальный язык программирования) можно быстро настроить управление объектами на сцене, также есть поддержка языка программирования Python.

 

 Рисунок 7 – Сцена в среде 3D движка Unreal Engine 4 с манипуляторами для тестирования кинематики
Рисунок 7 – Сцена в среде 3D движка Unreal Engine 4 с манипуляторами для тестирования кинематики

После проверки и отладки кода для расчета кинематики код был адаптирован для языка программирования C++ и перенесен в отдельную библиотеку, которая подгружается модулем управления роботом-манипулятором.

Рисунок 8 – Стенд системы с конвейером и роботом для захвата предметов
Рисунок 8 – Стенд системы с конвейером и роботом для захвата предметов

В ходе работы над проектом команде Zebrains удалось разработать и реализовать два собственных алгоритма обратной кинематики для робота xArm 2.0, адаптированного под задачу захвата объектов с конвейера с помощью вакуумной присоски. Первый алгоритм обеспечивает перемещение захвата в заданную точку при условии параллельности определённых звеньев, а второй — с учётом фиксированного угла наклона конечного звена. Оба решения реализованы на C++ и интегрированы в ROS2-ноды, что позволило добиться стабильного и предсказуемого поведения манипулятора в реальных условиях.

Полученные формулы и код открыты для обсуждения — особенно с учётом того, что подобные задачи часто возникают у других разработчиков, работающих с нестандартными или упрощёнными конфигурациями промышленных роботов. Мы надеемся, что описание наших подходов поможет вам либо напрямую использовать эти решения, либо адаптировать их под свои задачи. Если у вас есть идеи по улучшению, замечания по геометрической модели или опыт применения подобных методов на других манипуляторах - делитесь в комментариях.