【Matlab 六自由度机器人】求运动学逆解
往期回顾
【汇总】
相关资源：
【主线】
【补充说明】
- 关于灵活工作空间与可达工作空间的理解
- 关于改进型D-H参数(modified Denavit-Hartenberg)的详细建立步骤
- 关于旋转的参数化（欧拉角、姿态角、四元数）的相关问题
- 关于双变量函数atan2(x,y)的解释
- 关于机器人运动学反解的有关问题
前言
本篇介绍机器人逆运动学及逆解的求解方法，首先介绍如何理解逆向运动学，并利用D-H参数及正向运动学的齐次变换矩阵对机器人运动学逆解进行求解。
在博主读研期间，刚开始对机器人运动学逆解不甚了解，但在经历一次自己完全去钻研、尝试，独立写出属于自己的运动学逆解后，对机器人的结构有更加深刻的理解。因此希望初学机器人的朋友们也能够尝试独立去编写、尝试出自己的逆解代码，我相信，当你看到逆解解出的角度完全符合预期，也能够产生相当的成就感。
下面是运动学逆解的正文内容，主要讲述运动学逆解的方式和公式推算，最后进行代码的实现。
正文
一、运动学逆解
首先我们已经建立起了一个六自由度的机器人，在 D-H 参数 及机器人的 期望位姿 已知的情况下，求解机器人满足该期望位姿的各个关节变量，这种求解过程叫做机器人运动学逆解过程。求解机器人的逆运动学解法分为封闭解法和数值解法两种，其中封闭解法是指具有解析形式的解法，其求解速度比数值解法更快。
本文所研究的六自由度机器人，是具有六个旋转关节的串联型六自由度机器人，即RRR型的机器人。六自由度机器人的逆解封闭解法具体又可分为代数解法和几何解法，代数解法推导解析解的过程复杂，运算量较大，几何解法根据直观的几何关系获得关节变量的解析表达式，计算量较小。
本文主要讲解两种机器人运动学逆解的两种方法：
1. Pieper法
2. 常规法
1. Pieper 法
根据Pieper’s Solution的研究，如果六自由度机器人具有三个连续的轴交在同一点的情况，则手臂具有解析解。即六个关节均为旋转关节且后三个关节轴相交于一点的串联型机器人，其逆解具有封闭解。可以通过如下代数解法求解逆运动学的封闭解，该解法称为“三轴相交的 Pieper 解法”。下面进行各轴转角角度的公式推导。
我们知道，当后3个关节轴相交时，连杆坐标系 的原点位于该交点，该交点在基坐标系
中的坐标如下：
对于齐次变换矩阵的变换通式的理解请参阅【Matlab 六自由度机器人】运动学正解（附MATLAB机器人正解完整代码）
由改进D-H参数的连杆齐次变换矩阵式，应用齐次变换矩阵的变换通式
可知当 时，
应用齐次变换矩阵对于 可以得到
的表达式:
对 和
应用齐次变换矩阵得到
计算过程如下:
令 ，由上式可知
由上式写出 轴分量的方程，可得系统的两个方程:
对于，其计算如下:
之后求解 ，分三种情况:
(1) 若 ，则
，此时
已知。
的右边是关于
的函数。代入下式后，由包含
的一元二次方程可以解出
;
(2) 若 ，则
，此时
已知。代入上式后，利用上面的一元 二次方程解出
(3) 否则，从式中消去
和
，得到
代入解出
后，可得到一个一元四次方程，由此可解出
。 解出
后，再解出
，再根据
解出
。
Pieper 解法求解出了前三轴的关节变量，该解通常具有4组解或8组。由于后三轴的关节角只影响工具坐标系的方位，因此可以通过欧拉角解法求出后三轴的关节变量。
六自由度机器人工具坐标系的姿态，可以使用一组欧拉角 来表示，该转角序列表示绕
轴，
轴和
轴的旋转序列：先绕
轴旋转
角，再绕新坐标系的
轴旋转
角，最后绕新坐标系 的
轴旋转
角。欧拉角变换 Euler
可以通过三个旋转矩阵相乘表示
如果已知一个表示任意旋转的齐次变换, 即已知姿态矩阵
则根据姿态矩阵 可求出其对应的欧拉角:
对于6R型六自由度机器人，已经求得前三轴的关节变量，即可求出齐次变换矩 阵
假设建立欧拉角 Euler
的坐标系为
，则从基座标系
到坐标系
的齐次变换矩阵为:
式中，
一表示改进 DH参数表中的第四轴的连杆扭角。 欧拉角表示的资态矩阵的计算方法如式(2-21)所示:
其中
一表示齐次变换矩阵
中的旋转矩阵的转置；
一表示齐次变换矩阵
中的旋转矩阵。
求出欧拉角姿态矩阵后，可通过求出相应的欧拉角。所求得的欧 拉角与后三轴的关苖变量的关系为:
由于机器人腕关节翻转，可以得到后三轴的另外一组解:
上述算法可计算出的多组解，根据关节运动范围的限制可以舍弃一部分解，在余下的有效解中，通常选择最接近机器人当前位姿的“最短行程解”。
2. 《机器人学》常规求解
在求解运动学逆解前，我们先掌握运动学正解的知识，具体过程请参阅：【Matlab 六自由度机器人】运动学正解（附MATLAB机器人正解完整代码）
如前所述，矩阵右式的元素，或为零，或为常数，或者为第至第
个关节变量的函数。矩阵相等表明其对应元素分别相等，并可从每一矩阵方程得到
个方程式，每个方程式对应于
个矢量
的每一分量。我们先以之前建立过的 六自由度机器人 为例来阐述这些方程的求解。
该机器人的运动方程（设为方程式1）如下：
末端连杆的位姿已经给定，即
和
为已知，则求关节变量
的值称为运动反解。用末知的连杆逆变换左乘方程式1两边，把关节变量分离出来，从而求解。具体步骤如下:
1. 求
可用逆变换 左乘方程式1两边
令上式矩阵方程两瑞的元素 对应相等，可得
利用三角代换
式中， ，通过代换得到
的解。
式中，正、负号对应于 的两个可能解。
2. 求
在选定 的一个解之后，再令矩阵方程
两端的元素的
和
分别对应相等，即得两方程
式与
的平方和为
式中，
方程中已经消去
，且方程
与
具有相同形式，因而可由三角代换求解得到
式中，正、负号对应 的两种可能解。
3. 求
为求解 ，在矩阵方程1两边左乘逆变换
。
式中，变换 由正解出来的矩阵推出。令上式矩阵方程两边的元素
和 (2.4) 分别对应 相等可得
联立求解得
和
和
表达式的分母相等，且为正。于是
根据 和
解的四种可能组合，由上式可以得到相应的四种可能值
，于是可得到
的四种可能解
式中， 取与
相对应的值。
4. 求
因为式 (3.75) 的左边均为已知，令两边元素 和
分别对应相等，则可得
只要 ，便可求出
当 时，机械手处于奇㫒形位。此时，关节轴 4 和 6 重合，只能解出
与
的和或 差。奇异形位可以由式 (3.79)中
的两个变量是否都接近零来判别。若都接近零，则为奇异形位，否则，不是奇异形位。在奇异形位时，可任意选取
的值,再计箅相应的
值。
5. 求
据求出的 ，可进一步解出
，将方程1两端左乘逆变换
,
因下面方程式的左边 和
均已解出, 逆变换
为
上式方程式的右边 ，由正解得到的
给出。据矩阵两边元素
和
分别对应相等，可得
由此得到 的封闭解
6. 求
将方程式1改写为
令上式的矩阵方程两边元素 和
分别对应相等，可得
从而可求出 的封闭解
六自由度机器人的运动学反解可能存在
种解。但是，由于结构的限制，例如各关节变量不 能在全部
范围内运动，有些解不能实现。在机器人存在多种解的情况下，应选取其中满意的一组解，以满足机器人的工作要求。
后续更新简单理解的步骤
二、代码实现
1. 机器人模型的建立
根据往期文章，对六自由度机器人具体建模方法可参阅以下文章：
2. Pieper 法求六自由度机器人逆解
建立好机器人模型后，已知期望位姿的矩阵。将所求解的机器人的模型放入逆解脚本中，脚本的输入为期望位姿矩阵
，输出为角度
代码如下：
function [J] = MODikine(a)
%% MOD-DH参数
%定义连杆的D-H参数
%连杆偏移
d1 = 398;
d2 = -0.299;
d3 = 0;
d4 = 556.925;
d5 = 0;
d6 = 165;
%连杆长度
a1 = 0;
a2 = 168.3;
a3 = 650.979;
a4 = 156.240;
a5 = 0;
a6 = 0;
%连杆扭角
alpha1 = 0;
alpha2 = pi/2;
alpha3 = 0;
alpha4 = pi/2;
alpha5 = -pi/2;
alpha6 = pi/2;
%%
nx=a(1,1);ox=a(1,2);ax=a(1,3);px=a(1,4);
ny=a(2,1);oy=a(2,2);ay=a(2,3);py=a(2,4);
nz=a(3,1);oz=a(3,2);az=a(3,3);pz=a(3,4);
% 求解关节角1
theta1_1 = atan2(py,px)-atan2(d2,sqrt(px^2+py^2-d2^2));
theta1_2 = atan2(py,px)-atan2(d2,-sqrt(px^2+py^2-d2^2));
% 求解关节角3
m3_1 = (px^2+py^2+pz^2-a2^2-a3^2-d2^2-d4^2)/(2*a2);
theta3_1 = atan2(a3,d4)-atan2(m3_1,sqrt(a3^2+d4^2-m3_1^2));
theta3_2 = atan2(a3,d4)-atan2(m3_1,-sqrt(a3^2+d4^2-m3_1^2));
% 求解关节角2
ms2_1 = -((a3+a2*cos(theta3_1))*pz)+(cos(theta1_1)*px+sin(theta1_1)*py)*(a2*sin(theta3_1)-d4);
mc2_1 = (-d4+a2*sin(theta3_1))*pz+(cos(theta1_1)*px+sin(theta1_1)*py)*(a2*cos(theta3_1)+a3);
theta23_1 = atan2(ms2_1,mc2_1);
theta2_1 = theta23_1 - theta3_1;
ms2_2 = -((a3+a2*cos(theta3_1))*pz)+(cos(theta1_2)*px+sin(theta1_2)*py)*(a2*sin(theta3_1)-d4);
mc2_2 = (-d4+a2*sin(theta3_1))*pz+(cos(theta1_2)*px+sin(theta1_2)*py)*(a2*cos(theta3_1)+a3);
theta23_2 = atan2(ms2_2,mc2_2);
theta2_2 = theta23_2 - theta3_1;
ms2_3 = -((a3+a2*cos(theta3_2))*pz)+(cos(theta1_1)*px+sin(theta1_1)*py)*(a2*sin(theta3_2)-d4);
mc2_3 = (-d4+a2*sin(theta3_2))*pz+(cos(theta1_1)*px+sin(theta1_1)*py)*(a2*cos(theta3_2)+a3);
theta23_3 = atan2(ms2_3,mc2_3);
theta2_3 = theta23_3 - theta3_2;
ms2_4 = -((a3+a2*cos(theta3_2))*pz)+(cos(theta1_2)*px+sin(theta1_2)*py)*(a2*sin(theta3_2)-d4);
mc2_4 = (-d4+a2*sin(theta3_2))*pz+(cos(theta1_2)*px+sin(theta1_2)*py)*(a2*cos(theta3_2)+a3);
theta23_4 = atan2(ms2_4,mc2_4);
theta2_4 = theta23_4 - theta3_2;
% 求解关节角4
ms4_1=-ax*sin(theta1_1)+ay*cos(theta1_1);
mc4_1=-ax*cos(theta1_1)*cos(theta23_1)-ay*sin(theta1_1)*cos(theta23_1)+az*sin(theta23_1);
theta4_1=atan2(ms4_1,mc4_1);
ms4_2=-ax*sin(theta1_2)+ay*cos(theta1_2);
mc4_2=-ax*cos(theta1_2)*cos(theta23_2)-ay*sin(theta1_2)*cos(theta23_2)+az*sin(theta23_2);
theta4_2=atan2(ms4_2,mc4_2);
ms4_3=-ax*sin(theta1_1)+ay*cos(theta1_1);
mc4_3=-ax*cos(theta1_1)*cos(theta23_3)-ay*sin(theta1_1)*cos(theta23_3)+az*sin(theta23_3);
theta4_3=atan2(ms4_3,mc4_3);
ms4_4=-ax*sin(theta1_2)+ay*cos(theta1_2);
mc4_4=-ax*cos(theta1_2)*cos(theta23_4)-ay*sin(theta1_2)*cos(theta23_4)+az*sin(theta23_4);
theta4_4=atan2(ms4_4,mc4_4);
% 求解关节角5
ms5_1=-ax*(cos(theta1_1)*cos(theta23_1)*cos(theta4_1)+sin(theta1_1)*sin(theta4_1))-ay*(sin(theta1_1)*cos(theta23_1)*cos(theta4_1)-cos(theta1_1)*sin(theta4_1))+az*(sin(theta23_1)*cos(theta4_1));
mc5_1= ax*(-cos(theta1_1)*sin(theta23_1))+ay*(-sin(theta1_1)*sin(theta23_1))+az*(-cos(theta23_1));
theta5_1=atan2(ms5_1,mc5_1);
ms5_2=-ax*(cos(theta1_2)*cos(theta23_2)*cos(theta4_2)+sin(theta1_2)*sin(theta4_2))-ay*(sin(theta1_2)*cos(theta23_2)*cos(theta4_2)-cos(theta1_2)*sin(theta4_2))+az*(sin(theta23_2)*cos(theta4_2));
mc5_2= ax*(-cos(theta1_2)*sin(theta23_2))+ay*(-sin(theta1_2)*sin(theta23_2))+az*(-cos(theta23_2));
theta5_2=atan2(ms5_2,mc5_2);
ms5_3=-ax*(cos(theta1_1)*cos(theta23_3)*cos(theta4_3)+sin(theta1_1)*sin(theta4_3))-ay*(sin(theta1_1)*cos(theta23_3)*cos(theta4_3)-cos(theta1_1)*sin(theta4_3))+az*(sin(theta23_3)*cos(theta4_3));
mc5_3= ax*(-cos(theta1_1)*sin(theta23_3))+ay*(-sin(theta1_1)*sin(theta23_3))+az*(-cos(theta23_3));
theta5_3=atan2(ms5_3,mc5_3);
ms5_4=-ax*(cos(theta1_2)*cos(theta23_4)*cos(theta4_4)+sin(theta1_2)*sin(theta4_4))-ay*(sin(theta1_2)*cos(theta23_4)*cos(theta4_4)-cos(theta1_2)*sin(theta4_4))+az*(sin(theta23_4)*cos(theta4_4));
mc5_4= ax*(-cos(theta1_2)*sin(theta23_4))+ay*(-sin(theta1_2)*sin(theta23_4))+az*(-cos(theta23_4));
theta5_4=atan2(ms5_4,mc5_4);
% 求解关节角6
ms6_1=-nx*(cos(theta1_1)*cos(theta23_1)*sin(theta4_1)-sin(theta1_1)*cos(theta4_1))-ny*(sin(theta1_1)*cos(theta23_1)*sin(theta4_1)+cos(theta1_1)*cos(theta4_1))+nz*(sin(theta23_1)*sin(theta4_1));
mc6_1= nx*(cos(theta1_1)*cos(theta23_1)*cos(theta4_1)+sin(theta1_1)*sin(theta4_1))*cos(theta5_1)-nx*cos(theta1_1)*sin(theta23_1)*sin(theta4_1)+ny*(sin(theta1_1)*cos(theta23_1)*cos(theta4_1)+cos(theta1_1)*sin(theta4_1))*cos(theta5_1)-ny*sin(theta1_1)*sin(theta23_1)*sin(theta5_1)-nz*(sin(theta23_1)*cos(theta4_1)*cos(theta5_1)+cos(theta23_1)*sin(theta5_1));
theta6_1=atan2(ms6_1,mc6_1);
ms6_2=-nx*(cos(theta1_2)*cos(theta23_2)*sin(theta4_2)-sin(theta1_2)*cos(theta4_2))-ny*(sin(theta1_2)*cos(theta23_2)*sin(theta4_2)+cos(theta1_2)*cos(theta4_2))+nz*(sin(theta23_2)*sin(theta4_2));
mc6_2= nx*(cos(theta1_2)*cos(theta23_2)*cos(theta4_2)+sin(theta1_2)*sin(theta4_2))*cos(theta5_2)-nx*cos(theta1_2)*sin(theta23_2)*sin(theta4_2)+ny*(sin(theta1_2)*cos(theta23_2)*cos(theta4_2)+cos(theta1_2)*sin(theta4_2))*cos(theta5_2)-ny*sin(theta1_2)*sin(theta23_2)*sin(theta5_2)-nz*(sin(theta23_2)*cos(theta4_2)*cos(theta5_2)+cos(theta23_2)*sin(theta5_2));
theta6_2=atan2(ms6_2,mc6_2);
ms6_3=-nx*(cos(theta1_1)*cos(theta23_3)*sin(theta4_3)-sin(theta1_1)*cos(theta4_3))-ny*(sin(theta1_1)*cos(theta23_3)*sin(theta4_3)+cos(theta1_1)*cos(theta4_3))+nz*(sin(theta23_3)*sin(theta4_3));
mc6_3= nx*(cos(theta1_1)*cos(theta23_3)*cos(theta4_3)+sin(theta1_1)*sin(theta4_3))*cos(theta5_3)-nx*cos(theta1_1)*sin(theta23_3)*sin(theta4_3)+ny*(sin(theta1_1)*cos(theta23_3)*cos(theta4_3)+cos(theta1_1)*sin(theta4_3))*cos(theta5_3)-ny*sin(theta1_1)*sin(theta23_3)*sin(theta5_3)-nz*(sin(theta23_3)*cos(theta4_3)*cos(theta5_3)+cos(theta23_3)*sin(theta5_3));
theta6_3=atan2(ms6_3,mc6_3);
ms6_4=-nx*(cos(theta1_2)*cos(theta23_4)*sin(theta4_4)-sin(theta1_2)*cos(theta4_4))-ny*(sin(theta1_1)*cos(theta23_4)*sin(theta4_4)+cos(theta1_2)*cos(theta4_4))+nz*(sin(theta23_4)*sin(theta4_4));
mc6_4= nx*(cos(theta1_2)*cos(theta23_4)*cos(theta4_4)+sin(theta1_2)*sin(theta4_4))*cos(theta5_4)-nx*cos(theta1_2)*sin(theta23_4)*sin(theta4_4)+ny*(sin(theta1_2)*cos(theta23_4)*cos(theta4_4)+cos(theta1_2)*sin(theta4_4))*cos(theta5_1)-ny*sin(theta1_2)*sin(theta23_4)*sin(theta5_4)-nz*(sin(theta23_4)*cos(theta4_4)*cos(theta5_4)+cos(theta23_4)*sin(theta5_4));
theta6_4=atan2(ms6_4,mc6_4);
% 整理得到4组运动学非奇异逆解
theta_MOD1 = [
theta1_1,theta2_1,theta3_1,theta4_1,theta5_1,theta6_1;
theta1_2,theta2_2,theta3_1,theta4_2,theta5_2,theta6_2;
theta1_1,theta2_3,theta3_2,theta4_3,theta5_3,theta6_3;
theta1_2,theta2_4,theta3_2,theta4_4,theta5_4,theta6_4;
];
% 将操作关节‘翻转’可得到另外4组解
theta_MOD2 = ...
[
theta1_1,theta2_1,theta3_1,theta4_1+pi,-theta5_1,theta6_1+pi;
theta1_2,theta2_2,theta3_1,theta4_2+pi,-theta5_2,theta6_2+pi;
theta1_1,theta2_3,theta3_2,theta4_3+pi,-theta5_3,theta6_3+pi;
theta1_2,theta2_4,theta3_2,theta4_4+pi,-theta5_4,theta6_4+pi;
];
J = [theta_MOD1;theta_MOD2];
end
3. 常规法求六自由度机器人逆解
由于自写的代码过长，阅读较为困难，因此贴上部分逆解代码过程。
建立好机器人模型后，已知期望位姿的矩阵。将所求解的机器人的模型放入逆解脚本中，脚本的输入为期望位姿矩阵
，输出为角度
这里的代表的是机器人的模型，
是已知的期望位姿的矩阵，输出是各轴的角度
function [J] = Albert_MODikin(robot,a)
%% MOD-DH参数
%连杆偏移 %连杆长度 %连杆扭角
d1 = robot.d(1);a1 = robot.a(1);alpha1 = robot.alpha(1);
d2 = robot.d(2);a2 = robot.a(2);alpha2 = robot.alpha(2);
d3 = robot.d(3);a3 = robot.a(3);alpha3 = robot.alpha(3);
d4 = robot.d(4);a4 = robot.a(4);alpha4 = robot.alpha(4);
d5 = robot.d(5);a5 = robot.a(5);alpha5 = robot.alpha(5);
d6 = robot.d(6);a6 = robot.a(6);alpha6 = robot.alpha(6);
%%
nx=a(1,1);ox=a(1,2);ax=a(1,3);px=a(1,4);
ny=a(2,1);oy=a(2,2);ay=a(2,3);py=a(2,4);
nz=a(3,1);oz=a(3,2);az=a(3,3);pz=a(3,4);
%j1
j1a = py - ay*d6;
j1b = px - ax*d6;
j11 = atan2(j1a,j1b)-atan2(-d2, sqrt(j1a^2+j1b^2-d2^2));
j12 = atan2(j1a,j1b)-atan2(-d2,-sqrt(j1a^2+j1b^2-d2^2));
% 这部分是进行一个判断，将小于1e-16的数字看作0，建议删除。
if j11<1e-16
j11 = 0;
end
if j12<1e-16
j12 = 0;
end
%j3
j3a1 = px*cos(j11) - d6*(ax*cos(j11) + ay*sin(j11)) + py*sin(j11) - a2;
j3b1 = pz - d1 - az*d6;
j3k1 = j3a1^2 + j3b1^2 - a4^2 - d4^2 - a3^2;
j3d1 = j3k1/(2*a3);
j3a2 = px*cos(j12) - d6*(ax*cos(j12) + ay*sin(j12)) + py*sin(j12) - a2;
j3b2 = pz - d1 - az*d6;
j3k2 = j3a2^2 + j3b2^2 - a4^2 - d4^2 - a3^2;
j3d2 = j3k2/(2*a3);
j31 = atan2(j3d1, sqrt(abs(a4^2 + d4^2 - j3d1^2))) - atan2(a4,d4);%j11
j32 = atan2(j3d1,-sqrt(abs(a4^2 + d4^2 - j3d1^2))) - atan2(a4,d4);%j11
j33 = atan2(j3d2, sqrt(abs(a4^2 + d4^2 - j3d2^2))) - atan2(a4,d4);%j12
j34 = atan2(j3d2,-sqrt(abs(a4^2 + d4^2 - j3d2^2))) - atan2(a4,d4);%j12
% 经过一系列运算后，最终得到J
J = [j11 j21 j31 j41 j51 j61;
j11 j22 j32 j42 j52 j62;
j12 j23 j33 j43 j53 j63;
j12 j24 j34 j44 j54 j64;
j11 j21 j31 j45 j55 j65;
j11 j22 j32 j46 j56 j66;
j12 j23 j33 j47 j57 j67;
j12 j24 j34 j48 j58 j68];
end
程序运行后的结果如下：
总结
以上就是逆向运动学的内容，本文详细介绍了如何理解运动学逆解，以及对运动学逆解的求解方法的罗列及原理说明。最后对运动学逆解实现代码的编程，使逆解能够有效得到目标位姿所对应的各轴转角。
参考文献
文章出处登录后可见！
