机械臂正向与逆向运动学求解_机械臂求逆解

        机械臂的正运动学求解即建立DH参数表，然后计算出各变换矩阵以及最终的变换矩阵。逆运动学求解，即求出机械臂各关节θ角与px,py,pz的关系，建立θ角与末端姿态之间的数学模型，在这里以IRB6700为例，对IRB6700进行正逆运动学求解和验证。

目录

正运动学求解

逆运动学求解

正逆运动学模型的验证

正运动学验证

逆运动学验证

总的Matlab代码，包含正逆运动学求解和验证

参考文献

---

正运动学求解

        首先使用DH法建立坐标系如下：

         查阅IRB6700的参数如下表

根据变换矩阵公式编写matlab代码计算得到正运动学的各变换矩阵

上式中的为,为，为

 计算变换矩阵的MATLAB代码如下

clear;clc
%%
%导入参数
%注，这里的a,afa的下标在实际使用时均需要减去1
syms a afa d theta [1 6]
a = [0 a2 a3 a4 0 0];
afa = [0 -90 0 -90 90 -90];
d = [d1 0 0 d4 0 d6];
theta = [theta1 theta2 theta3 theta4 theta5 theta6];
syms T [4 4 6]
%%
%正运动学模型求解
%计算各变换矩阵及总的变换矩阵
for i = 1:6
   T(:,:,i) =  trans_cal(  afa(i), a(i), d(i), theta(i)*180/pi );
          if(i == 1)
              trans_matrix = T(:,:,i);
          else
              trans_matrix = trans_matrix*T(:,:,i);
          end
end

 用到的funcion函数如下

function T = trans_cal(afa_ii,a_ii,d_i,theta_i)
%%
%计算变换矩阵函数T_{i-1,i}
%输入的参数为，afa_{i-1},a_{i-1},d_i,theta_i,与DH表达法的参数表对应
%ii 为i-1
%注意，这里输入的角度，均采用角度制，不采用弧度制
 
T = [cosd(theta_i)              -sind(theta_i)             0               a_ii
     sind(theta_i)*cosd(afa_ii)  cosd(theta_i)*cosd(afa_ii)  -sind(afa_ii)    -sind(afa_ii)*d_i
    sind(theta_i)*sind(afa_ii)   cosd(theta_i)*sind(afa_ii)  cosd(afa_ii)     cosd(afa_ii)*d_i
     0                         0                         0                1               ];
 
end

 各个变换矩阵存储在T中，其中即为变换矩阵,存储在trans_matrix中。

逆运动学求解

设总的变换矩阵设为下数：

有等式

根据上面的正运动学求解可以得到矩阵中的等各元素的值分别为什么

使用matlab的实时计算脚本可以较具体地看到

 然后通过封闭解法求逆，将上述等式中右边地变换矩阵乘到左边，然后观察等式两边矩阵中地元素，使等式两边的矩阵中的某个元素相等，不断分离变量，然后求解得到角

 首先求解，构造等式

 

 观察等式，使等式两边矩阵的第三行第四列元素相等

上述等式只含有一个未知数，一个等式一个未知数，求解过程如下：

注：为了便于计算，${\color{Red}a _{i-1} }$均用${\color{Red}a _{i} }$代替，${\color{Red}\alpha_{i-1} }$也用${\color{Red}\alpha_{i} }$代替，即下列式子中的${\color{Red}a _{1} }$实际上为${\color{Red}a _{0} }$，其他的${\color{Red}a }$和${\color{Red}\alpha}$亦然如此

 用同样的方式求解$\theta_{2}$

 令等式左右两边的第一行第四列，和第三行第四列分别相等，求解过程如下

 

然后是$\theta_{3}$

 令等式两边第一行第四列，第二行第四列相等

求解过程如下：

 再然后是$\theta_{4}$的求解

令左右两边第三行第三列相等

求解过程如下

 最后是$\theta_{5}$和$\theta_{6}$的求解

对$\theta_{5}$ 构造，使等式左右第一行第三列，和第三行第三列相等

对$\theta_{6}$构造，使等式左右第三行第一列，和第三行第二列相等

 求解过程如下

 至此完成了所有$\theta$角的求解

 再次注名：上述推导过程中，为了便于计算，${\color{Red}a _{i-1} }$均用${\color{Red}a _{i} }$代替，${\color{Red}\alpha_{i-1} }$也用${\color{Red}\alpha_{i} }$代替，即上述式子中的${\color{Red}a _{1} }$实际上为${\color{Red}a _{0} }$，其他的${\color{Red}a }$和${\color{Red}\alpha}$亦然如此

正逆运动学模型的验证

正运动学验证

随机选定一组关节角$\theta_{1} \sim \theta_{6}$进行验证，

使用之前求解的正运动学模型求解得到末端姿态矩阵$_{6}^{0}\textrm{T}$如下

使用RobotStudio软件导入IRB6700机器人，输入选择的$\theta_{1} \sim \theta_{6}$角

注意：由连杆参数可知，RobotStudio软件中的$\theta_{2}$角度为原角度加上$\frac{\pi }{2}$，$\theta_{6}$为原角度减去$\pi$，所以RobotStudio中输入的角度应如下：

 

得到的TCP末端xyz坐标       

 

对比RobotStudio软件中的末端姿态[x,y,z]值与正解求得的矩阵中末端姿态左边[px,py,pz]可得，两者相等，因此正运动学模型正确。

逆运动学验证

由上述的正运动学验证继续计算，上述正运动学得到的末端姿态矩阵如下：

根据此末端姿态矩阵对逆运动学进行求解，可得到八种不同的$\theta$角组合，使末端姿态为[1635.672,594.531,1396.902]

求解的matlab代码较长，共一个主m文件，以及5个function文件，放在本文最后。

求解得到共八种组合的$\theta$角结果，根据IRB6700机器人的各关节角工作范围减去4个不符合的解后，得到四个符合的$\theta$角组合，然后保持RobotStudio软件中的$\theta$角参数不变，调整RobotStudio软件中的轴配置参数，将逆运动学模型求解得到的结果与RobotStudio软件中的结果进行对比如下表

由上表的结果对比可知，除了计算精度问题导致的，部分结果小数点后两位与RobotStudio存在一个单位的差异外，结果一致，因为可以验证逆运动学模型正确。

总的Matlab代码，包含正逆运动学求解和验证

clear;clc
%%
%导入参数
%注，这里的a,afa的下标在实际使用时均需要减去1
syms a afa d theta [1 6]
a = [0 a2 a3 a4 0 0];
afa = [0 -90 0 -90 90 -90];
d = [d1 0 0 d4 0 d6];
theta = [theta1 theta2 theta3 theta4 theta5 theta6];
syms T [4 4 6]
%%
%正运动学模型求解
%计算各变换矩阵及总的变换矩阵
for i = 1:6
   T(:,:,i) =  trans_cal(  afa(i), a(i), d(i), theta(i)*180/pi );
          if(i == 1)
              trans_matrix = T(:,:,i);
          else
              trans_matrix = trans_matrix*T(:,:,i);
          end
end
 
 
%%
%逆运动学计算求解
syms r [3 3]
syms PX_X PY_Y PZ_Z
T_60 = [r1_1 r1_2 r1_3 PX_X;
        r2_1 r2_2 r2_3 PY_Y;
        r3_1 r3_2 r3_3 PZ_Z;
        0    0    0    1];
    
left_1 = simplify( inv( T(:,:,2) ) * inv( T(:,:,1) ) * T_60 * inv( T(:,:,6) ) );
right_1 = simplify( T(:,:,3)*T(:,:,4)*T(:,:,5) );
 
 
left_2 =  simplify( inv( T(:,:,1) )* T_60 * inv( T(:,:,6) ) );
right_2 = simplify( T(:,:,2)*T(:,:,3)*T(:,:,4)*T(:,:,5) );
 
left_3 =  simplify( inv(   T(:,:,1)*T(:,:,2)   ) * T_60 * inv( T(:,:,6) ) );
right_3 = simplify( T(:,:,3)*T(:,:,4)*T(:,:,5) );
 
left_4 =  simplify(  inv(   T(:,:,1)*T(:,:,2)*T(:,:,3)*T(:,:,4)*T(:,:,5)   ) * T_60 );
right_4 = simplify(  T(:,:,6)  );
 
left_5 =  simplify(  inv(   T(:,:,1)*T(:,:,2)*T(:,:,3)*T(:,:,4) ) * T_60 );
right_5 = simplify(  T(:,:,5)*T(:,:,6)  );
 
left_6 =  simplify(  inv(   T(:,:,1)*T(:,:,2)*T(:,:,3)*T(:,:,4)*T(:,:,5)   ) * T_60 );
right_6 = simplify(  T(:,:,6)  );
 
 
%%
%逆运动学验算
 
%末端姿态
T_ni = [-0.5    0   0.866   1635.672;
        0       1   0       594.531;
        -0.866  0   -0.5    1396.902;
        0       0   0       1];
    
    
r1_1 = T_ni(1,1);r1_2 = T_ni(1,2);r1_3 = T_ni(1,3);PX_X   = T_ni(1,4);
r2_1 = T_ni(2,1);r2_2 = T_ni(2,2);r2_3 = T_ni(2,3);PY_Y   = T_ni(2,4);
r3_1 = T_ni(3,1);r3_2 = T_ni(3,2);r3_3 = T_ni(3,3);PZ_Z   = T_ni(3,4);
 
syms a d [1 6]
d1 = 780;d4 = 1142.5;d6 = 200;
a2 = 320;a3 = 1125  ;a4 = 200; 
 
syms A B C D E F [1 6]
 
%theta1的解
theta1_1 = ( atan2(0,1) - atan2( d6*r2_3-PY_Y , PX_X-d6*r1_3 ) )*180/pi %22.1230
theta1_2 = ( -atan2(0,-1) - atan2( d6*r2_3-PY_Y , PX_X-d6*r1_3 ) )*180/pi %-157.8770
 
 
 
%theta2的解
[theta2_1 , theta2_2] = theta2_calculate(theta1_1)  %theta1为22.123°时
[theta2_3 , theta2_4] = theta2_calculate(theta1_2)  %theta1为-157.8770°时
 
 
%theta3的解
theta3_11_21    =   theta3_calculate(theta1_1,theta2_1)  %theta1为22.123°theta2为-81.3955°时
theta3_11_22    =   theta3_calculate(theta1_1,theta2_2)  %theta1为22.123°theta2为22.0674°时
theta3_12_23    =   theta3_calculate(theta1_2,theta2_3)  %theta1为-157.8770°theta2为172.8834°时
theta3_12_24    =   theta3_calculate(theta1_2,theta2_4)  %theta1为-157.8770°theta2为228.4872°时
 
 
 
%theta4的解
[theta4_1_11_21,theta4_2_11_21] = theta4_calculate(theta1_1,theta2_1,theta3_11_21)  %theta1为22.123°theta2为-81.3955°时
[theta4_1_11_22,theta4_2_11_22] = theta4_calculate(theta1_1,theta2_2,theta3_11_22)  %theta1为22.123°theta2为22.0674°时
[theta4_1_12_23,theta4_2_12_23] = theta4_calculate(theta1_2,theta2_3,theta3_12_23)  %theta1为-157.8770°theta2为172.8834°时
[theta4_1_12_24,theta4_2_12_24] = theta4_calculate(theta1_2,theta2_4,theta3_12_24)  %theta1为-157.8770°theta2为228.4872°时
 
 
%theta5的解
theta5_11_21_41 = theta5_calculate(theta1_1,theta2_1,theta3_11_21,theta4_1_11_21) %theta1为22.123°theta2为-81.3955°,theta4为第1种
theta5_11_21_42 = theta5_calculate(theta1_1,theta2_1,theta3_11_21,theta4_2_11_21) %theta1为22.123°theta2为-81.3955°,theta4为第2种
theta5_11_22_41 = theta5_calculate(theta1_1,theta2_2,theta3_11_22,theta4_1_11_22) %theta1为22.123°theta2为22.0674°,theta4为第1种
theta5_11_22_42 = theta5_calculate(theta1_1,theta2_2,theta3_11_22,theta4_2_11_22) %theta1为22.123°theta2为22.0674°,theta4为第2种
 
theta5_12_23_41 = theta5_calculate(theta1_2,theta2_3,theta3_12_23,theta4_1_12_23) %theta1为-157.8770°theta2为172.8834°,theta4为第1种
theta5_12_23_42 = theta5_calculate(theta1_2,theta2_3,theta3_12_23,theta4_2_12_23) %theta1为-157.8770°theta2为172.8834°,theta4为第2种
theta5_12_24_41 = theta5_calculate(theta1_2,theta2_4,theta3_12_24,theta4_1_12_24) %theta1为-157.8770°theta2为228.4872°,theta4为第1种
theta5_12_24_42 = theta5_calculate(theta1_2,theta2_4,theta3_12_24,theta4_2_12_24) %theta1为-157.8770°theta2为228.4872°,theta4为第2种
 
 
% theta6的解
theta6_11_21_41 = theta6_calculate(theta1_1,theta2_1,theta3_11_21,theta4_1_11_21) %theta1为22.123°theta2为-81.3955°,theta4为第1种
theta6_11_21_42 = theta6_calculate(theta1_1,theta2_1,theta3_11_21,theta4_2_11_21) %theta1为22.123°theta2为-81.3955°,theta4为第2种
theta6_11_22_41 = theta6_calculate(theta1_1,theta2_2,theta3_11_22,theta4_1_11_22) %theta1为22.123°theta2为22.0674°,theta4为第1种
theta6_11_22_42 = theta6_calculate(theta1_1,theta2_2,theta3_11_22,theta4_2_11_22) %theta1为22.123°theta2为22.0674°,theta4为第2种
theta6_12_23_41 = theta6_calculate(theta1_2,theta2_3,theta3_12_23,theta4_1_12_23) %theta1为-157.8770°theta2为172.8834°,theta4为第1种
theta6_12_23_42 = theta6_calculate(theta1_2,theta2_3,theta3_12_23,theta4_2_12_23) %theta1为-157.8770°theta2为172.8834°,theta4为第2种
theta6_12_24_41 = theta6_calculate(theta1_2,theta2_4,theta3_12_24,theta4_1_12_24) %theta1为-157.8770°theta2为228.4872°,theta4为第1种
theta6_12_24_42 = theta6_calculate(theta1_2,theta2_4,theta3_12_24,theta4_2_12_24) %theta1为-157.8770°theta2为228.4872°,theta4为第2种
 

function [theta2_1,theta2_2] = theta2_calculate(theta1)
%%
%导入参数
    syms r [3,3]
    syms a afa d
    syms PX_X PY_Y_Y PZ_Z
    T_ni = [-0.5    0   0.866   1635.672;
        0       1   0       594.531;
        -0.866  0   -0.5    1396.902;
        0       0   0       1];    
    r1_1 = T_ni(1,1);r1_2 = T_ni(1,2);r1_3 = T_ni(1,3);PX_X   = T_ni(1,4);
    r2_1 = T_ni(2,1);r2_2 = T_ni(2,2);r2_3 = T_ni(2,3);PY_Y   = T_ni(2,4);
    r3_1 = T_ni(3,1);r3_2 = T_ni(3,2);r3_3 = T_ni(3,3);PZ_Z   = T_ni(3,4);
    syms a d [1 6]
    d1 = 780;d4 = 1142.5;d6 = 200;
    a2 = 320;a3 = 1125  ;a4 = 200; 
    syms A B C D E F [1 6]
    theta1 = theta1/180*pi;
 
%%
%theta2计算
    A2 = ( PX_X - d6*r1_3 )*cos(theta1) + ( PY_Y - d6*r2_3 )*sin(theta1) - a2;
    B2 = d1 + d6*r3_3 - PZ_Z;
    C2 = 2*A2*a3;
    D2 = 2*B2*a3;
    E2 = A2^2 + B2^2 + a3^2 - a4^2 - d4^2;
    F2 = sqrt( C2^2 + D2^2 );
    
    theta2_1 = ( -atan2(C2,D2)  + atan2(  E2/F2 ,  sqrt( 1- (E2/F2)^2 ) ) )*180/pi  ;
    theta2_2 = ( -atan2(C2,D2)  + atan2(  E2/F2 , -sqrt( 1- (E2/F2)^2 ) ) )*180/pi  ;
 
end

function theta3 = theta3_calculate(theta1,theta2)
%%
%导入参数
    syms r [3,3]
    syms a afa d
    syms PX_X PY_Y_Y PZ_Z
    T_ni = [-0.5    0   0.866   1635.672;
        0       1   0       594.531;
        -0.866  0   -0.5    1396.902;
        0       0   0       1];    
    r1_1 = T_ni(1,1);r1_2 = T_ni(1,2);r1_3 = T_ni(1,3);PX_X   = T_ni(1,4);
    r2_1 = T_ni(2,1);r2_2 = T_ni(2,2);r2_3 = T_ni(2,3);PY_Y   = T_ni(2,4);
    r3_1 = T_ni(3,1);r3_2 = T_ni(3,2);r3_3 = T_ni(3,3);PZ_Z   = T_ni(3,4);
    syms a d [1 6]
    d1 = 780;d4 = 1142.5;d6 = 200;
    a2 = 320;a3 = 1125  ;a4 = 200; 
    syms A B C D E F [1 6]
    theta1 = theta1/180*pi;
    theta2 = theta2/180*pi;
 
%%
%theta3计算
    A3 = d1*sin(theta2) - a2*cos(theta2) - PZ_Z*sin(theta2) + d6*r3_3*sin(theta2) + PX_X*cos(theta1)*cos(theta2) ...
    + PY_Y*cos(theta2)*sin(theta1) - d6*r1_3*cos(theta1)*cos(theta2) - d6*r2_3*cos(theta2)*sin(theta1) - a3;
 
    B3 = d1*cos(theta2) - PZ_Z*cos(theta2) + a2*sin(theta2) + d6*r3_3*cos(theta2) - PX_X*cos(theta1)*sin(theta2) ...
    - PY_Y*sin(theta1)*sin(theta2) + d6*r2_3*sin(theta1)*sin(theta2) + d6*r1_3*cos(theta1)*sin(theta2);
 
    C3 = ( a4*B3-d4*A3 ) / ( a4^2 + d4^2 );
    
    D3 = ( a4*A3+d4*B3 ) / ( a4^2 + d4^2 );
    
    theta3 = atan2( C3,D3 )*180/pi;
 
 
end

function [theta4_1,theta4_2] = theta4_calculate(theta1,theta2,theta3)
%%
%导入参数
    syms r [3,3]
    syms a afa d
    syms PX_X PY_Y_Y PZ_Z
    T_ni = [-0.5    0   0.866   1635.672;
        0       1   0       594.531;
        -0.866  0   -0.5    1396.902;
        0       0   0       1];    
    r1_1 = T_ni(1,1);r1_2 = T_ni(1,2);r1_3 = T_ni(1,3);PX_X   = T_ni(1,4);
    r2_1 = T_ni(2,1);r2_2 = T_ni(2,2);r2_3 = T_ni(2,3);PY_Y   = T_ni(2,4);
    r3_1 = T_ni(3,1);r3_2 = T_ni(3,2);r3_3 = T_ni(3,3);PZ_Z   = T_ni(3,4);
    syms a d [1 6]
    d1 = 780;d4 = 1142.5;d6 = 200;
    a2 = 320;a3 = 1125  ;a4 = 200; 
    syms A B C D E F [1 6]
    theta1 = theta1/180*pi;
    theta2 = theta2/180*pi;
    theta3 = theta3/180*pi;
 
%%
%theta4计算
 
A4 = r3_3*cos(theta2)*sin(theta3) + r3_3*cos(theta3)*sin(theta2) - r1_3*cos(theta1)*cos(theta2)*cos(theta3) ...
    - r2_3*cos(theta2)*cos(theta3)*sin(theta1) + r1_3*cos(theta1)*sin(theta2)*sin(theta3) ...
    + r2_3*sin(theta1)* sin(theta2)*sin(theta3);  
B4 = r2_3*cos(theta1) - r1_3*sin(theta1);
theta4_1 = ( atan2( 0,1 ) + atan2( B4,A4 ) )*180/pi;
theta4_2 = ( atan2( 0,-1 ) + atan2( B4,A4 ) )*180/pi;
 
end

function theta5 = theta5_calculate(theta1,theta2,theta3,theta4)
%%
%导入参数
    syms r [3,3]
    syms a afa d
    syms PX_X PY_Y_Y PZ_Z
    T_ni = [-0.5    0   0.866   1635.672;
        0       1   0       594.531;
        -0.866  0   -0.5    1396.902;
        0       0   0       1];    
    r1_1 = T_ni(1,1);r1_2 = T_ni(1,2);r1_3 = T_ni(1,3);PX_X   = T_ni(1,4);
    r2_1 = T_ni(2,1);r2_2 = T_ni(2,2);r2_3 = T_ni(2,3);PY_Y   = T_ni(2,4);
    r3_1 = T_ni(3,1);r3_2 = T_ni(3,2);r3_3 = T_ni(3,3);PZ_Z   = T_ni(3,4);
    syms a d [1 6]
    d1 = 780;d4 = 1142.5;d6 = 200;
    a2 = 320;a3 = 1125  ;a4 = 200; 
    syms A B C D E F [1 6]
    theta1 = theta1/180*pi;
    theta2 = theta2/180*pi;
    theta3 = theta3/180*pi;
    theta4 = theta4/180*pi;
 
%%
%theta4计算
 
A5 = r1_3*sin(theta1)*sin(theta4) - r2_3*cos(theta1)*sin(theta4) - r3_3*cos(theta2)*cos(theta4)*sin(theta3) ...
    - r3_3*cos(theta3)*cos(theta4)*sin(theta2) + r1_3*cos(theta1)*cos(theta2)*cos(theta3)*cos(theta4) + ...
    r2_3*cos(theta2)*cos(theta3)*cos(theta4)*sin(theta1) - r1_3*cos(theta1)*cos(theta4)*sin(theta2)*sin(theta3) ...
    - r2_3*cos(theta4)*sin(theta1)*sin(theta2)*sin(theta3);
B5 = r3_3*sin(theta2)*sin(theta3) - r3_3*cos(theta2)*cos(theta3) - r1_3*cos(theta1)*cos(theta2)*sin(theta3) ... 
    - r1_3*cos(theta1)*cos(theta3)*sin(theta2) - r2_3*cos(theta2)*sin(theta1)*sin(theta3) - ... 
    r2_3*cos(theta3)*sin(theta1)*sin(theta2);
theta5 = atan2(-A5,B5)*180/pi;
 
end

function theta6 = theta6_calculate(theta1,theta2,theta3,theta4)
%%
%导入参数
    syms r [3,3]
    syms a afa d
    syms PX_X PY_Y_Y PZ_Z
    T_ni = [-0.5    0   0.866   1635.672;
        0       1   0       594.531;
        -0.866  0   -0.5    1396.902;
        0       0   0       1];    
    r1_1 = T_ni(1,1);r1_2 = T_ni(1,2);r1_3 = T_ni(1,3);PX_X   = T_ni(1,4);
    r2_1 = T_ni(2,1);r2_2 = T_ni(2,2);r2_3 = T_ni(2,3);PY_Y   = T_ni(2,4);
    r3_1 = T_ni(3,1);r3_2 = T_ni(3,2);r3_3 = T_ni(3,3);PZ_Z   = T_ni(3,4);
    syms a d [1 6]
    d1 = 780;d4 = 1142.5;d6 = 200;
    a2 = 320;a3 = 1125  ;a4 = 200; 
    syms A B C D E F [1 6]
    theta1 = theta1/180*pi;
    theta2 = theta2/180*pi;
    theta3 = theta3/180*pi;
    theta4 = theta4/180*pi;
 
%%
%theta4计算
 
A6 = r2_1*cos(theta1)*cos(theta4) - r1_1*cos(theta4)*sin(theta1) - r3_1*cos(theta2)*sin(theta3)*sin(theta4) ...
    - r3_1*cos(theta3)*sin(theta2)*sin(theta4) + r1_1*cos(theta1)*cos(theta2)*cos(theta3)*sin(theta4) + ...
    r2_1*cos(theta2)*cos(theta3)*sin(theta1)*sin(theta4) - r1_1*cos(theta1)*sin(theta2)*sin(theta3)*sin(theta4) ...
    - r2_1*sin(theta1)*sin(theta2)*sin(theta3)*sin(theta4);
B6 = r2_2*cos(theta1)*cos(theta4) - r1_2*cos(theta4)*sin(theta1) - r3_2*cos(theta2)*sin(theta3)*sin(theta4) ... 
    - r3_2*cos(theta3)*sin(theta2)*sin(theta4) + r1_2*cos(theta1)*cos(theta2)*cos(theta3)*sin(theta4)  ...
    + r2_2*cos(theta2)*cos(theta3)*sin(theta1)*sin(theta4) - r1_2*cos(theta1)*sin(theta2)*sin(theta3)*sin(theta4) ...
    - r2_2*sin(theta1)*sin(theta2)*sin(theta3)*sin(theta4);
 
theta6 = atan2(-A6,-B6)*180/pi;
 
end

function T = trans_cal(afa_ii,a_ii,d_i,theta_i)
%%
%计算变换矩阵函数T_{i-1,i}
%输入的参数为，afa_{i-1},a_{i-1},d_i,theta_i,与DH表达法的参数表对应
%ii 为i-1
%注意，这里输入的角度，均采用角度制，不采用弧度制
 
T = [cosd(theta_i)              -sind(theta_i)             0               a_ii
     sind(theta_i)*cosd(afa_ii)  cosd(theta_i)*cosd(afa_ii)  -sind(afa_ii)    -sind(afa_ii)*d_i
    sind(theta_i)*sind(afa_ii)   cosd(theta_i)*sind(afa_ii)  cosd(afa_ii)     cosd(afa_ii)*d_i
     0                         0                         0                1               ];
 
end

参考文献

[1] 尹绪伟. 打磨机器人不同位姿下的刚度特性研究[D]. 武汉: 武汉理工大学, 2019.

[2] 王宇迪. 基于加工姿态最优的发动机飞轮壳机器人铣削修边空间路径规划[D]. 武汉: 武汉理工大学, 2022.