在推导机器人运动学逆解的解析解时,经常会遇到以下三角方程:
k 1 s i n ( θ ) + k 2 c o s ( θ ) = k 3 (1) k_1sin(\theta)+k_2cos(\theta)=k_3 \tag{1} k1sin(θ)+k2cos(θ)=k3(1)
求解以上三角方程的解析解对运动学逆解的推导过程至关重要。下面采用两种方法进行推导。
一、推导步骤
1、方法1
式(1)联立以下三角恒等式:
s
i
n
2
(
θ
)
+
c
o
s
2
(
θ
)
=
1
(2)
sin^2(\theta)+cos^2(\theta)=1 \tag{2}
sin2(θ)+cos2(θ)=1(2)
利用MATLAB解符号方程组(1)(2),代码如下:
clc;
clear;
syms k1 k2 k3 theta sinTheta cosTheta real
result = solve(k1 * sinTheta + k2 * cosTheta == k3, sinTheta^2 + cosTheta^2 == 1, [sinTheta, cosTheta]);
sinTheta = simplify(result.sinTheta)
cosTheta = simplify(result.cosTheta)
很容易求得两组解:
{
s
i
n
θ
=
(
k
1
k
3
∓
k
2
k
1
2
+
k
2
2
−
k
3
2
)
/
(
k
1
2
+
k
2
2
)
c
o
s
θ
=
(
k
2
k
3
±
k
1
k
1
2
+
k
2
2
−
k
3
2
)
/
(
k
1
2
+
k
2
2
)
(3)
\begin{cases} sin\theta= (k_1k_3 \mp k_2\sqrt{k_1^2 + k_2^2 - k_3^2})/(k_1^2 + k_2^2) \\ cos\theta=(k_2k_3 \pm k_1\sqrt{k_1^2 + k_2^2 - k_3^2})/(k_1^2 + k_2^2) \\ \tag 3 \end{cases}
{sinθ=(k1k3∓k2k12+k22−k32)/(k12+k22)cosθ=(k2k3±k1k12+k22−k32)/(k12+k22)(3)
若
k
1
2
+
k
2
2
<
k
3
2
k_1^2+k_2^2<k_3^2
k12+k22<k32,则当前给定位姿区域不可到达;若
k
1
=
k
2
=
0
k_1=k_2=0
k1=k2=0,则当前给定位姿为奇异位姿。这里,采用代数的方法便可方便判断给定位姿是否在机器人可达工作空间里,并判断该位姿是否为奇异位姿!
若
k
1
2
+
k
2
2
≥
k
3
2
k_1^2+k_2^2\ge k_3^2
k12+k22≥k32,由于
k
1
2
+
k
2
2
>
0
k_1^2+k_2^2>0
k12+k22>0,易得:
θ
=
a
t
a
n
2
(
k
1
k
3
∓
k
2
k
1
2
+
k
2
2
−
k
3
2
,
k
2
k
3
±
k
1
k
1
2
+
k
2
2
−
k
3
2
)
(4)
\theta=atan2(k_1k_3 \mp k_2\sqrt{k_1^2 + k_2^2 - k_3^2},k_2k_3 \pm k_1\sqrt{k_1^2 + k_2^2 - k_3^2}) \tag{4}
θ=atan2(k1k3∓k2k12+k22−k32,k2k3±k1k12+k22−k32)(4)
这里采用双变量反正切函数
a
t
a
n
2
atan2
atan2,有几个优点,详见博文:为什么机器人运动学逆解最好采用双变量反正切函数atan2而不用反正/余弦函数?
2、方法2
令
t
a
n
(
θ
/
2
)
=
u
tan(\theta/2)=u
tan(θ/2)=u,将三角恒等式
c
o
s
θ
=
(
1
−
u
2
)
/
(
1
+
u
2
)
,
s
i
n
θ
=
2
u
/
(
1
+
u
2
)
cos\theta=(1-u^2)/(1+u^2),sin\theta=2u/(1+u^2)
cosθ=(1−u2)/(1+u2),sinθ=2u/(1+u2)代入式(1),整理得:
(
k
2
+
k
3
)
u
2
−
2
k
1
u
+
k
3
−
k
2
=
0
(5)
(k_2+k_3)u^2-2k_1u+k_3-k_2=0 \tag{5}
(k2+k3)u2−2k1u+k3−k2=0(5)
当
k
2
+
k
3
=
0
k_2+k_3=0
k2+k3=0时,
θ
=
±
π
\theta=\pm \pi
θ=±π。
当
k
2
+
k
3
≠
0
k_2+k_3\ne0
k2+k3=0时,若
k
1
2
+
k
2
2
<
k
3
2
k_1^2+k_2^2<k_3^2
k12+k22<k32,则当前给定位姿区域不可到达;
若
k
1
2
+
k
2
2
≥
k
3
2
k_1^2+k_2^2\ge k_3^2
k12+k22≥k32,求解式(5)关于
u
u
u的一元二次方程的两实根,得到三角方程(1)的解:
θ
=
2
a
t
a
n
(
(
k
1
±
k
1
2
+
k
2
2
−
k
3
2
)
/
(
k
2
+
k
3
)
)
(6)
\theta= 2atan((k_1 \pm \sqrt{k_1^2+k_2^2-k_3^2})/(k_2+k_3))\tag{6}
θ=2atan((k1±k12+k22−k32)/(k2+k3))(6)
二、实例(以SCARA机器人逆解为例)
在博文:scara机器人运动学正逆解中,采用几何法推导得到机器人的运动学逆解,本文采用代数法来求运动学逆解。
SCARA机器人的运动学正解为:
{
x
=
L
1
c
o
s
θ
1
+
L
2
c
o
s
(
θ
1
+
θ
2
)
y
=
L
1
s
i
n
θ
1
+
L
2
s
i
n
(
θ
1
+
θ
2
)
(7)
\left \{ \begin{array}{c} x=L_1cos\theta_1+L_2cos(\theta_1+\theta_2) \\ \tag 7 y=L_1sin\theta_1+L_2sin(\theta_1+\theta_2) \end{array}\right.
{x=L1cosθ1+L2cos(θ1+θ2)y=L1sinθ1+L2sin(θ1+θ2)(7)
式(7)可消去
s
i
n
(
θ
1
+
θ
2
)
,
c
o
s
(
θ
1
+
θ
2
)
sin(\theta_1+\theta_2),cos(\theta_1+\theta_2)
sin(θ1+θ2),cos(θ1+θ2),转化为关于
θ
1
\theta_1
θ1的三角方程:
2
y
L
1
s
i
n
θ
1
+
2
x
L
1
c
o
s
θ
1
=
x
2
+
y
2
+
L
1
2
−
L
2
2
(8)
2yL_1sin\theta_1+2xL_1cos\theta_1=x^2+y^2+ L_1^2-L_2^2\tag{8}
2yL1sinθ1+2xL1cosθ1=x2+y2+L12−L22(8)
令
k
1
=
2
y
L
1
,
k
2
=
2
x
L
1
,
k
3
=
x
2
+
y
2
+
L
1
2
−
L
2
2
k_1=2yL_1,k_2=2xL_1,k_3=x^2+y^2+ L_1^2-L_2^2
k1=2yL1,k2=2xL1,k3=x2+y2+L12−L22,可解得:
θ
1
=
a
t
a
n
2
(
k
1
k
3
∓
k
2
k
1
2
+
k
2
2
−
k
3
2
,
k
2
k
3
±
k
1
k
1
2
+
k
2
2
−
k
3
2
)
(9)
\theta_1=atan2(k_1k_3 \mp k_2\sqrt{k_1^2 + k_2^2 - k_3^2},k_2k_3 \pm k_1\sqrt{k_1^2 + k_2^2 - k_3^2}) \tag{9}
θ1=atan2(k1k3∓k2k12+k22−k32,k2k3±k1k12+k22−k32)(9)
由式(7)得:
{
s
i
n
(
θ
1
+
θ
2
)
=
(
y
−
L
1
s
i
n
θ
1
)
/
L
2
c
o
s
(
θ
1
+
θ
2
)
=
(
x
−
L
1
c
o
s
θ
1
)
/
L
2
(10)
\left \{ \begin{array}{c} sin(\theta_1+\theta_2)=(y-L_1sin\theta_1)/L_2 \\ \tag {10} cos(\theta_1+\theta_2)=(x-L_1cos\theta_1)/L_2 \end{array}\right.
{sin(θ1+θ2)=(y−L1sinθ1)/L2cos(θ1+θ2)=(x−L1cosθ1)/L2(10)
解得:
θ
2
=
a
t
a
n
2
(
y
−
L
1
s
i
n
θ
1
,
x
−
L
1
c
o
s
θ
1
)
−
θ
1
(11)
\theta_2=atan2(y-L_1sin\theta_1,x-L_1cos\theta_1)-\theta_1\tag{11}
θ2=atan2(y−L1sinθ1,x−L1cosθ1)−θ1(11)
MATLAB代码:
clc;
clear;
L1 = 200.0; %mm
L2 = 300.0; %mm
theta1 = 1.0; %rad
theta2 = 3.0; %rad
x = L1 * cos(theta1) + L2 * cos(theta1 + theta2);
y = L1 * sin(theta1) + L2 * sin(theta1 + theta2);
disp(['输入角度(rad):',num2str(theta1), ',', num2str(theta2)])
disp(['输入x,y(mm):',num2str(x), ',', num2str(y)])
k1 = 2.0 * y * L1;
k2 = 2.0 * x * L1;
k3 = x^2 + y^2 + L1^2 - L2^2;
temp = k1^2 + k2^2 - k3^2;
if temp < -eps
disp('区域不可到达');
return;
elseif temp < eps
disp('奇异点');
else
%do nothing
end
theta1 = zeros(2, 1);
theta1(1) = atan2(k1 * k3 - k2 * sqrt(temp), k2 * k3 + k1 * sqrt(temp));
theta1(2) = atan2(k1 * k3 + k2 * sqrt(temp), k2 * k3 - k1 * sqrt(temp));
s1 = sin(theta1);
c1 = cos(theta1);
theta2 = atan2(y - L1 * s1, x - L1 * c1) - theta1;
for i = 1 : 2
if theta2(i) > pi
theta2(i) = theta2(i) - 2.0 * pi;
end
if theta2(i) < -pi
theta2(i) = theta2(i) + 2.0 * pi;
end
end
xx = L1 * cos(theta1) + L2 * cos(theta1 + theta2);
yy = L1 * sin(theta1) + L2 * sin(theta1 + theta2);
disp(['输出角度(rad):',num2str(theta1(1)), ',', num2str(theta2(1)), '; ', num2str(theta1(2)), ',', num2str(theta2(2))])
disp(['输出x,y(mm):',num2str(xx(1)), ',', num2str(yy(1)), '; ', num2str(xx(2)), ',', num2str(yy(2))])