Yummy Robot Development Notes
机械臂模型(MDH)
机械臂设计符合Pieper判据,即相邻的三个关节轴相交于一点,根据改进DH法可得:
i |
θi |
di(m) |
ai−1(m) |
αi−1(rad) |
1 |
θ1 |
0(0.11) |
0 |
0 |
2 |
θ2 |
0 |
0 |
α1=π/2 |
3 |
θ3 |
0 |
$a_2=$0.3 |
0 |
4 |
θ4 |
0.27 |
$a_3=$0.096 |
α3=π/2 |
5 |
θ5 |
0 |
0 |
α4=−π/2 |
6 |
θ6 |
0(0.107) |
0 |
α5=π/2 |
工作空间
使用MATLAB,对机械臂蒙特卡洛采样法可建立笛卡尔空间坐标系下的可视范围:
正运动学
正运动学就是通过各个关节的角度和连杆间的变换关系解算出末端的位姿
0T6=⎣⎢⎢⎢⎡nxnynz0oxoyoz0axayaz0pxpypz1⎦⎥⎥⎥⎤=0T1(θ1)1T2(θ2)2T3(θ3)3T4(θ4)4T5(θ5)5T6(θ6)
相邻两关节之间的变换:
ii−1T=⎣⎢⎢⎢⎡cθisθicαi−1sθisαi−10−sθicθicαi−1cθisαi−100−sαi−1cαi−10ai−1−sαi−1dicαi−1di1⎦⎥⎥⎥⎤
逆运动学
逆解分析
机械臂逆解的个数
串联机械臂的逆运动学问题求解相比正运动学复杂度更高,对于6自由度机械臂来说,在逆运动学问题分析过程中需要对齐次变换矩阵中的相互独立的6个非线性超越方程进行求解。该方程有6个未知量,其解的存在性以及解的个数由机械臂的运动学结构、目标点是否位于灵活工作空间以及程序中队关节限制而决定的。
在存在可行解的情况下,其个数不唯一,对于6自由度机械臂而言,最多可能存在16个解,解的个数与DH参数中的a(连杆长度)有关:
ai |
逆解个数 |
a1=a3=a5=0 |
≤4 |
a3=a5=0 |
≤8 |
a3=0 |
≤16 |
All ai=0 |
≤16 |
对于满足pieper判据的6自由度机器人来说,其后三轴相交于一点,可看作一个球型的腕关节,其前三个关节和后三个关节可以相互解耦,由下图可以看出前三轴决定了腕关节点空间中的位置,而后三个关节则是空间中的旋转方向,前三关节有以下4种情况,而整个腕关节一般有正反两种情况(J4,J6反转互补),因此一共有8种解。

机械臂逆解类型
机械臂逆运动学分为封闭解法和数值解法,封闭解法通过解析方式求封闭解,分为代数法和几何法得到的解称为解析解;数值法通过迭代求解,速度较慢
Pieper判据
判断机械臂是否存在解析解
在求解机械臂逆运动学前判断其解的存在性非常重要。现有研究表明,所有的串联6自由度机械臂均是可解的,但是这种解只能通过数值解法得到,计算难度大,复杂度高。因此在确定逆运动学解法前,需要探究机械臂逆运动学问题是否存在解析解,Pieper判据是机器人领域常用的一种判断解析解是否存在的判据。
Pipper判据:6自由度机械臂存在解析解的充分条件是相邻的三个关节旋转轴相交于一点
机械臂逆解求解
简易机械臂通用求解
几何法:
逆运动学就是根据设定的笛卡尔空间位姿反求关节空间中的各个关节的角度,如上图所示,我们设定目标位置(x,y)和最后一共关节的朝向角度ϕ,通过三角函数求解θ1,θ2,θ3。
我们首先连接两点构建一个三角形,可以用余弦定理求解θ2(橙色):
x2+y2=l12+l22−2l1l2cos(180−θ2)
我们便可以表示出:
cos(θ2)=2l1l2x2+y2−l12−l22
我们再利用另一个内角ψ,使用余弦定理:
cos(ψ)=−2l1x2+y2l22−(x2+y2)−l12
由于ψ是三角形内角所以满足0<ψ<180,我们可以通过ψ求解θ1,但这个时候出现了如上图所示多解的情况:
θ1={β+ψβ−ψθ2<0θ2>0
θ1={atan2(y,x)+ψatan2(y,x)−ψθ2<0θ2>0
(相对角度逆时针为正,顺时针为负)
最后设定角度ϕ为三个角度加和:
θ1+θ2+θ3=ϕ
便可求解出θ3:
θ3=ϕ−θ1−θ2
代数解法:
根据连杆参数很容易求得这个机械臂的运动学方程:
TWB=T30⎣⎢⎢⎢⎡c123s12300−s123c123000010l1c1+l2c12l1s1+l2s1201⎦⎥⎥⎥⎤
其中:
c123=cos(θ1+θ2+θ3)
同几何法一样,我们设定目标位置(x,y)和ϕ,同样可以得出变换矩阵:
WBT=⎣⎢⎢⎢⎡cϕsϕ00−sϕcϕ000010xy01⎦⎥⎥⎥⎤
所有的可达点都要位于上式的子空间中,令上下两式子相等,可以求得4个非线性方程:
c_{\phi}=& c_{123} &&\text{(1)}\\
s_{\phi}=& s_{123} &&\text{(2)}\\
x =& l_1c_1+l_2c_{12} &&\text{(3)}\\
y =& l_1s_1+l_2s_{12} &&\text{(4)}\\
(3)、(4)同时平方然后相加得:
x2+y2=l12+l22+2l1l2c2
便可求解出c2:
c2=2l1l2x2+y2−l12−l22
上述有解的条件是右边的值必须在[−1,1]。这个约束可以用来检查解是否存在,若约束条件不满足,则目标点位置太原,操作臂不可达。
假定目标点在工作空间中,s2表达式为:
s2=±1−c22
最后应用双变量反正切公式计算θ2,得:
θ2=Atan2(s2,c2)
上式的正负号选择应用于多解,可以选择“肘部朝上”或“肘部朝下”,在确定θ2时,应用循环方法求解运动学参数,先确定期望关节角的正弦和余弦,然后应用双变量反正切公式的方法。这样我们确保得出所有的解,且所求的角度是在适当的象限里。
求出了θ2可以根据(3)(4)式求出θ1,可以将(3)(4)写成以下形式:
x=k1c1−k2s1y=k1s1+k2c1
其中:
\begin{alignat}{2}
&k_1=&&l_1+l_2c_2 \\
&k_2=&&l_2s_2
\end{alignat}
为了求解这种形式的方程,可进行变量代换,实际上就是改变常数k1和k2的形式
令
\begin{alignat}{2}
&r=+\sqrt{k_1^2+k_2^2} \\
&\gamma=Atan2(k_2,k_1)
\end{alignat}
则:
k1=rcosγk2=rsinγ
因此可改写为:
rx=cosγcosθ1−sinγsinθ1ry=cosγsinθ1+sinγcosθ1
进而可得:
cos(γ+theta1)=rxsin(γ+theta1)=ry
利用反正切公式可得:
γ+θ1=Atan2(ry,rx)=Atan2(y,x)
从而
θ1=Atan2(y,x)−Atan2(k2,k1)
注意,θ2符号的选取也将导致k2符号的变化,因此影响到θ1。利用万能公式进行变化求解的方法经常出现在运动学求解过程中。
最后可由θ1,θ2,θ3的和:
θ1+θ2+θ3=Atan2(sϕ,cϕ)=ϕ
得出θ3:
θ3=ϕ−θ1−θ2
三轴相交的Pieper解法
尽管一般的6自由度机器人没有封闭解,但是在某些特殊情况下还是可解的。本项目中机械臂也是按照Pieper判据设计,即三个连续轴相交于一点的6自由度操作臂,(包括具有3个连续平行轴的机械臂,可以认为交点在无穷远处)Pieper也提出了针对6个旋转关节且后面3个轴相交的操作臂的逆运动学求解方法。
当最后的三个轴相交时,连杆坐标系{4}、{5}、{6}的原点均位于这个交点上:
这点在基坐标系中的位置是:
0P4ORG=10T21T32T3P4ORG=⎣⎢⎢⎢⎡xyz1⎦⎥⎥⎥⎤
在MDH下,我们可以用下式表示相邻两个关节间的变换关系:
ii−1T=⎣⎢⎢⎢⎡cθisθicαi−1sθisαi−10−sθicθicαi−1cθisαi−100−sαi−1cαi−10ai−1−sαi−1dicαi−1di1⎦⎥⎥⎥⎤
当i=4时,第四列即为3P4ORG,则上式可化为:
0P4ORG=10T21T32T⎣⎢⎢⎢⎡a3−d4sα3d4cα31⎦⎥⎥⎥⎤
在ii−1T中,除了θi是需要求解的量,其余全都可以是DH参数已经确定的,因此每个ii−1T都可看作关于θi的函数:
0P4ORG=10T21T⎣⎢⎢⎢⎡f1(θ3)f2(θ3)f3(θ3)1⎦⎥⎥⎥⎤
其中:
⎣⎢⎢⎢⎡f1(θ3)f2(θ3)f3(θ3)1⎦⎥⎥⎥⎤=32T⎣⎢⎢⎢⎡a3−d4sα3d4cα31⎦⎥⎥⎥⎤
展开得出f1,f2,f3:
\begin{align}
\MoveEqLeft
f_1=a_3c_3 + d_4s\alpha_3s_3+a_2\\
f_2=&a_3c\alpha_2s_3 - d_4s\alpha_3c\alpha_2c_3 -d_4s\alpha_2c\alpha_3-d_3s\alpha_2 \\
f_3=&a_3s\alpha_2s_3-d_4s\alpha_3s\alpha_2c_3+d_4c\alpha_2c\alpha_3+d_3c\alpha_2
\end{align}
看起来复杂,但是代入本项目中的参数后为:
\begin{alignat}{2}
&f_1=&&0.096*c_3+0.27s_3+0.3 \\
&f_2=&&0.096*s_3-0.27*c_3-0.27 \\
&f_3=&&0
\end{alignat}
同样地,我们对10T,21T应用变换得:
0P4ORG=⎣⎢⎢⎢⎡c1g1−s1g2s1g1+c1g2g31⎦⎥⎥⎥⎤(1)
其中:
\begin{alignat}{2}
&g_1=c_2f_1-s_2f_2+a_1 \\
&g_2=s_2c\alpha_1f_1+c_2c\alpha_1f_2-s\alpha_1f_3-d_2s\alpha_1 \\
&g_3=s_2s\alpha_1f_1+c_2s\alpha_1f_2-c\alpha_1f_3-d_2c\alpha_1 \\
\end{alignat}
\begin{alignat}{2}
&g_1=c_2f_1-s_2f_2 \\
&g_2=0 \\
&g_3=s_2f_1+c_2f_2 \\
\end{alignat}
其中值得注意的是10T只有旋转矩阵部分,平移向量为0,因为按照MDH的坐标系建立方法两者坐标系在同一个点上。
现在写出0P4ORG绝对值平方的表达式,这里r=x2+y2+z2:
r=g12+g22+g32
展开得:
r=f12+f22+f32+a12+d22+2d2f3+2a1(c2f1−s2f2)
r=f12+f22
现在,写出Z方向的分量,可化为以下两个方程:
r=(k1c2+k2s2)2a1+k3z=(k1s2−k2c2)sα1+k4(2)
其中:
\begin{alignat}{2}
&k_1=f_1 \\
&k_2=-f_2 \\
&k_3=f_1^2+f_2^2+f_3^2+a_1^2+d_2^2+2d_2f_3 \\
&k_4=f_3c\alpha_1+d_2c\alpha_1
\end{alignat}
上式消去了因变量θ1,并且因变量θ2的关系式简单。
现讨论如何求解θ3,分三种情况:
- 若a1=0,则r=k3,这里r是已知的。右边(k3)仅是关于θ3的函数。代入下式:
\begin{alignat}{2}
u=&tan\frac{\theta}{2} \\
cos\theta=&\frac{1-u^2}{1+u^2} \tag{*}
\end{alignat}
由包含tan2θ3的二次方程可以解出θ3。
- 若sα1=0,则z=k4,这里z是已知的,再次代入(*)式后,利用上面的一元二次方程可以解出θ3。
- 否则,从式(2)中消去s2和c2,得到:
4a12(r−k3)2+s2α1(z−k4)2=k12+k22
再次代入(*)式,可以得到一个4次方程,由此可解出θ3。
末端角度解算
为了完成求解工作,还需要求出θ4,θ5,θ6。由于这些轴相交,故这些关节角只影响末端连杆的方向,我们只需要60R的旋转分量就能计算出这三个角度。在求出前三个关节角度后,可以由θ4=0时连杆坐标系{4}相对于基坐标系的方向计算出04R∣θ4=0。坐标系{6}的期望方向于连杆坐标系{4}的方向的差别仅在于最后三个关节的作用。由于60R已知,因此这个问题可以通过以下计算得出结果:
64R∣θ4=0=40R−1∣θ4=060R
对于大多数机械臂来说,都能够按照Z−Y−Z欧拉角解法解出最后三个关节角,对于任何一个4、5、6轴相交的机械臂来说,最后三个关节角能够通过一组合适的欧拉角来定义。
按照ZYZ的旋转顺序,分别定义:
α=θ4β=θ5γ=θ6
可得到其旋转矩阵:
64RZYZ(α,β,γ)=⎣⎢⎡cαcβcγ−sαsγsαcβcγ+cαsγ−sβcγ−cαcβsγ−sαcγ−sαcβsγ+cαcγsβsγcαsβsαsβcβ⎦⎥⎤
令:
64RZYZ(α,β,γ)=⎣⎢⎡r11r21r31r12r22r32r13r23r33⎦⎥⎤
如果sinβ=0,可得到:
\begin{alignat}{2}
\beta=&Atan2(\sqrt{r_{31}^2+r_{32}^2},r_{33}) \\
\alpha=&Atan2(r_{23}/s\beta,r_{13}/s\beta) \\
\gamma=&Atan2(r_{32}/s\beta,-r_{31}/s\beta)
\end{alignat}
上式只是取了β的正平方根,还可以有第二个解,但是我们总是求满足β∈[0,π]的单解。
如果sinβ=0,即θ5=0或者θ5=π,这时J4和J6共线,上式的解就退化了,在这种情况下,仅能求出α和β的和或者差,在这种情况下一般取α=0,结果如下:
如果β=0,则解为:
\begin{alignat}{2}
\beta=&0 \\
\alpha=&0 \\
\gamma=&Atan2(-r_{12},r_{11})
\end{alignat}
如果β=π,则解为:
\begin{alignat}{2}
\beta=&180 \\
\alpha=&0 \\
\gamma=&Atan2(r_{12},-r_{11})
\end{alignat}
雅可比矩阵
雅可比矩阵推导
雅可比矩阵是多维形式的导数。例如,假设有6个函数,每个函数都有6个独立的变量:
y1=f1(x1,x2,x3,x4,x5,x6)y2=f2(x1,x2,x3,x4,x5,x6)...y6=f6(x1,x2,x3,x4,x5,x6)
也可以用矢量符号表示这些等式:
Y=F(X)
要计算出yi的微分关于xi的微分的函数,可简单应用多元函数求导法则计算:
δY=δXδFδX
上式中的6X6偏微分矩阵就是雅可比矩阵。如果f1(X)到f6(x)都是非线性函数,那么这些偏微分都是xi的函数,因此可以表达为:
δY=J(X)δX
同时将上式两端除以时间的微分,可以将雅可比矩阵看成X中的速度向Y中速度的映射:
Y˙=J(X)X˙
在任一时刻,X都有一个确定的值,J(X)是一个线性变换。在每一个新的时刻,如果X变化,线性变换也会随之改变。所以,雅可比是时变的线性变换。
在机器人学中,通常使用雅可比矩阵将关节速度于操作比末端的笛卡尔速度联系起来:
0v=0J(Θ)Θ˙
上式中,Θ是机械臂关节角矢量,v是笛卡尔速度矢量。上式中左上标代表笛卡尔速度所参考的坐标系。对于任意已知的机械臂位形,关节速度和机械臂末端速度的关系是线性的,然而这种线性关系仅是瞬间的,因为在下一时刻,雅可比矩阵就会有微小的变化。
对于6关节机器人,雅可比矩阵是6X6维的,Θ是6X1维的,0v是6X1维的,由线速度矢量和角速度矢量排列起来的:
0v=[0v0ω]
雅可比矩阵的行数等于机械臂在笛卡尔空间中自由度的数量,列数等于机械臂的关节数量。
奇异性
已知一个线性变换可以将关节速度和笛卡尔速度联系起来,那么这个线性变换矩阵是否可逆呢?这个矩阵是否是非奇异的呢?如果是非奇异的,那么如果已知笛卡尔空间中的速度,就可以对该矩阵求逆计算出关节空间的速度:
Θ˙=J−1(Θ)v
但是有一个问题就是,由于雅可比矩阵随着关节空间位置不同也随之改变,是否对于所有Θ都是可逆的?如果不是,在什么位置不可逆?
大多数机械臂都有使得雅可比矩阵奇异的Θ值,这些位置就是奇异点。所有的机械臂柱子工作空间的边界都存在奇异点,而且在它们工作空间内也有奇异点,大致分为两类:
-
工作空间边界的奇异点,出现在机械臂完全展开或收回使得末端执行器处于或非常接近工作空间边界的情况。
-
工作空间内部的奇异点,总是远离工作空间的边界,通常由于两个或两个以上的关节轴线共线引起的。
当机械臂处于奇异点时,它会失去一个或多个自由度(笛卡尔空间),这也就是说,在笛卡尔空间的某个方向上(或某个子空间中),无论选择什么样的关节速度,都不能使得机器人手臂运动。
利用雅可比矩阵检测奇异点
-
雅可比矩阵的秩,在奇异点时,雅可比矩阵的秩降低,导致自由度丢失。
- 奇异点个数:min(6,n)−rank(J)
-
奇异点检测:
计算行列式det(J),越接近0,则接近奇异点
-
当det(J)=0,Jb=不可逆(秩亏缺),机械臂处于奇异点,某些方向运动能力丧失
-
当det(J)≈0,J接近奇异点,关节速度将会趋近于无穷大,或数值不稳定,不可控。
-
行列式代表了矩阵变换的"体积缩放因子",行列式为零意味着雅可比矩阵将高维空间压缩到了低维空间。
可以通过MATLAB将某点的雅可比矩阵JJT绘制成椭球,椭球越“圆”,表示机械臂在各方向上运动能力越均衡。 椭球越“瘪”或“压扁”,说明在某些方向运动能力差。下图分别为线速度和角速度的可操作性椭球。
或计算条件数cond(J),越大越接近奇异点
-
条件数κ(J)衡量矩阵数值的稳定性:
κ(J)=σminσmax
其中σmax和σmin是J的最大和最小奇异值,可通过对J奇异值分解(SVD)求出当前时刻的奇异值。
-
当 $\kappa(J) \to \infty $即 $ \sigma_{\text{min}} \to 0 ,矩阵接近奇异,实际应用中,当\kappa(J)>10^3$时认为接近奇异。
-
条件数的优点:行列式受矩阵尺度的影响,单位变化导致行列式值的变化,但是条件数无量纲更稳定,可以直接反应矩阵求逆的误差放大倍数。
当在在某些姿态时,就可以根据条件数的曲线值来判断这条路径上产生奇异点的位置,如下图就是在沿X轴运动过程中出现的腕部奇异点附近的曲线:
奇异点规避
奇异点规避可以从设计阶段就增加冗余的自由度绕过奇异点,也可以对关节进行限位约束,避免进入奇异点,但是减少了部分工作空间。
也可以通过**阻尼最小二乘法(DLS)**改进:
-
求J的伪逆J+,通过奇异值分解(SVD):
J=UΣVT⇒J+=VΣ+UT
其中:
Σ=⎝⎜⎜⎜⎜⎜⎛σ1σ2...σr0⎠⎟⎟⎟⎟⎟⎞
Σ+=⎝⎜⎜⎜⎜⎜⎛σ11σ21...σr10⎠⎟⎟⎟⎟⎟⎞
-
此时可以放弃阈值小于某个值的奇异值,但会产生不连续性,利用阻尼最小二乘法计算伪逆的方式是在倒置奇异值之前为其设置一个下界:
J∗=JT(JJT+λI)−1
通过引入阻尼参数λ,避免分母为0,来进行平滑过渡:
J∗=V⎝⎜⎜⎜⎜⎜⎜⎛σ12+λσ10σ22+λσ20...0σr2+λσr0⎠⎟⎟⎟⎟⎟⎟⎞UT
最后得到:
Θ˙=JT(JJT+λI)−1v
效果就是λ越大,关节速度越平滑,但跟踪误差增加,适合实时控制。