引言

随着国内市场需求的日益提升和机器人技术的不断发展，机械臂已逐渐替代人力成为自动生产线上的主导力量^[1]，如在喷涂、焊接、零件装配、机械加工等作业中^[2-4].然而，机械臂的发展也面临许多挑战，其中之一在于机械臂对于负载能力的要求逐渐提高，如家具搬运^[5]，水果采摘^[6]等，而大负载的出现会导致许多问题.例如，大负载的出现会对机械臂的运动性能带来影响，大惯性力使机械臂的负载效应变得明显，从而影响其末端的定位精度^[7，8].因此，研究大负载如何影响机械臂末端精度成为一个极为重要的科学问题.

机械臂种类繁多，若按机械臂轴数来分，常见的有四轴、六轴、七轴.其中，四轴机械臂具有出色的刚性和高速度，但灵活性相对较低，适用于高速和高重复性的任务；六轴机械臂则具有高度的灵活性，但需要较大的运动空间，且随着连杆数量的增加，其动力学方程也变得更加复杂，适用于喷涂、焊接等任务；与六轴机械臂相较之下，七轴机械臂在灵活性和精度方面表现更加出色，同时运动空间小，但其刚性较低，逆运动学问题求解也更为复杂，因此，七轴机械臂适用于空间有限但需要高度灵活性的任务^[9-11].

机械臂的特性分析主要包括运动学分析和动力学分析.运动学分析涉及正、逆运动学方程的求解.其中正运动学主要采用DH（Denavit-Hartenberg）参数法得到连杆变换矩阵从而得到正运动学方程^[12，13]，此外，还可以通过螺旋理论^[14]、Gibbs-Appell公式^[15]建立运动学方程.而逆运动学常用的方法可分为解析法和数值法，对于复杂的机械臂一般找不到其解析解只能采用数值法，常用的数值法有雅可比迭代法^[16]，粒子群法^[17]等.对于动力学分析，常用的动力学建模方法有Lagrange法^[18]、Kane法^[19]、Newton-Euler法^[20]等.

机械臂的轨迹规划主要分为两种，一种是在关节空间中，一种是在笛卡尔空间中^[21].其中关节空间的轨迹规划计算方便简单，但其末端运动路径不直观，一般用于对末端执行器运动路径需求小的作业，方法有多项式插值法^[22]，RRT算法（Rapidly Exploring Random Tree）^[23]；而笛卡尔空间的轨迹规划计算困难复杂，但能直观地看到末端运动路径，一般用于对末端执行器运动路径有需求的作业，方法有遗传算法^[24]，梯形加减速法^[25]等.

要保证大负载下末端的移动性能则需要对机械臂的控制进行补偿.对于大负载所引起的关节变形的控制补偿，毛晨涛等^[26]通过激光跟踪仪测量辨识大负载机器人的关节刚度系数，并对控制律进行优化设计.王旭等^[27]提出了通过直接计算关节变形量补偿变形误差的方法.

对于大负载引起的惯性力的增加，罗天洪等^[8]提出了基于神经网络补偿的机械臂轨迹控制策略.李佳兵等^[28]利用PSO（Particle Swarm Optimization）优化算法对大负载机械臂载荷参数自动控制进行设计.对于大负载所引起的残余振动，马驰聘、张春林等^[29]通过施加速度负反馈控制的方法，有效地抑制系统的残余振动.Liu等^[30]提出了一种关节空间多目标轨迹规划方案来减小承载有效载荷柔性机械臂的残余振动.Xiao^[31]针对变负载提出了一种H_∞自适应控制策略来保证位置跟踪.

本文以七轴串联机械臂为对象，开展运动学分析、路径规划、动力学分析、稳定性分析、定位控制误差分析等研究工作，旨在探索大负载对机械臂末端精度的影响.针对这一研究，本文主要开展以下几项工作：（1）机械臂运动学建模及路径规划；（2）.含大负载的机械臂动力学建模；（3）控制律的设计与稳定性分析；（4）控制参数和末端负载对机械臂末端移动性能的影响.

1 七轴串联机械臂的运动学分析

机械臂的动力学分析是对机械臂系统施加精确控制的基础，而运动学分析则是建立动力学模型和施加动力学控制的前提条件.因此，本节将进行Franka七轴串联机械臂的运动学分析.

1.1 M-DH参数法

DH参数法是串联机械臂中最常用的一种运动学建模方法，目前，主要的DH参数法有两种，一个是标准的DH参数法，一个是改进的DH参数法（M-DH法（Modified Denavit-Hartenberg）），其二者的区别主要在于四个参数的选取标准不同.本文采用的建模方法为M-DH参数法.

如图1所建立的M-DH坐标系，其中，连杆坐标系{i-1}的z_i-1轴与i-1关节的轴线重合，方向任意；x_i-1轴与z_i-1轴与z_i轴之间的公垂线重合，若z_i-1轴与z_i轴相交，则x_i-1轴垂直于z_i-1轴与z_i轴所形成的平面；y_i-1轴由右手法则确定；原点o_i-1建立在x_i-1轴和z_i-1轴的交点，若z_i-1轴和z_i轴相交，则原点建立在z_i-1轴与z_i轴的交点处.

对于基底坐标系（{0}坐标系），其与机械臂的底端固定，方向可任意选择.为了方便起见，一般与当关节1的旋转角度为0时的连杆坐标系{1}重合.

对于末端坐标系（{n}坐标系），其原点和x_n轴的方向可以任意选择，为了方便起见，对于原点的选择一般使连杆参数d_n=0，对于x_n轴的选择一般使连杆参数θ_n=0.

由此可以得到相邻连杆间的坐标变换矩阵为：

其中参数a_i-1为{i-1}坐标系的z_i-1轴沿着x_i-1轴到{i}坐标系z_i轴的距离，参数α_i-1为{i-1}坐标系的z_i-1轴绕着x_i-1轴转到{i}坐标系z_i轴的角度，参数d_i为{i-1}坐标系的x_i-1轴沿着z_i轴到{i}坐标系x_i轴的距离，参数θ_i为{i-1}坐标系的x_i-1轴绕着z_i轴转到{i}坐标系x_i轴的角度.

1.2 正运动学分析

本文以Franka七轴串联机械臂为对象，其模型如图2所示.

根据Franka七轴串联机械臂模型图，建立其M-DH坐标系如图3所示.

下表是该机械臂的DH参数：

其中d₁=0.333m，d₃=0.315m，d₅=0.383m，d₇=0.159m，a₄=0.083m，a₅=-0.083m，a₇=-0.089m.

将DH参数代入式（1）中，求出相邻连杆间的变换矩阵后，依次相乘可得到Franka七轴串联机械臂的正运动学方程，如下：

1.3 逆运动学分析

机械臂的逆运动学分析是根据机械臂末端相对于机械臂基底坐标系的位姿来求解出机械臂每个关节的旋转角度.因此，本小节将采用雅可比迭代法，对Franka七轴串联机械臂的逆运动学方程进行求解.

根据机械臂的末端笛卡尔速度和关节速度微分运动方程，有：

其中J表示机械臂的雅可比矩阵，而dX，dq的形式如下：

则可得到：

将机械臂末端的微小运动dX用各关节的运动dq叠加来近似得到.dX越小，其近似效果越好，即认为机械臂各关节的旋转角度由初始角度q一步一步旋转到期望角度q_e，假设每一步的步长为Δq，有：

因为七轴串联机械臂属于冗余机械臂，因此雅可比矩阵求逆需引入伪逆矩阵，即：

其中pinv（*）表示求伪逆.

由此，求逆运动学解的具体步骤为：

a.设置好最大迭代次数，初始末端位姿矩阵和期望末端位姿矩阵，设定好位姿误差范数g；

b.计算下一步机械臂的各关节变换角度；

c.更新机械臂各关节角度q，其中q=q+Δq；

d.更新此时机械臂末端位姿状态X;

e.计算此时误差dX的1范数（‖X‖）；

f.迭代次数+1；

g.判断‖dX‖≤g，若是，则退出迭代，若不是，则返回步骤b，直至‖dX‖≤g或达到最大迭代次数；

h.输出此时各关节的角度q，q即为期望末端位姿的逆运动学解.

1.4 关节空间轨迹规划

本小节采用三次多项式插值法对机械臂关节空间进行轨迹规划，其函数θ（t）的形式如下：

由约束条件：

其中，起始位置的关节角度θ₀是已知的，而期望位置的关节角度θ_d可通过给定末端位姿后，通过逆运动学方程求解得到.而起始位置的关节角速度和期望位置的关节角速度由自己给定，也是已知的.因此，可得到：

2 七轴串联机械臂的动力学分析

2.1 含负载Franka七轴串联机械臂动力学建模

机械臂包含多个连杆和关节，其动力学系统具有多输入多输出的性质，存在复杂的耦合关系和严重的非线性.因此，本节以含负载的Franka七轴串联机械臂为对象，采用拉格朗日法对其进行动力学分析.

对于任意机械系统，规定其拉格朗日函数的形式为系统的动能L与势能U之差，即：

则该系统的动力学方程可写成：

其中动能（采用极惯性矩方法求解），势能分别为：

极惯性矩阵表示为：

其中，I_o=I_xx+I_yy+I_zz表示连杆相对原点的惯性矩.重力矢量g=[0  0  9.8  0]m/s².

推导出机械臂系统的总动能和总势能后，将其代回式（12）可得：

其中D（q）是机械臂的质量矩阵，C（q）是机械臂的科氏矩阵，G（q）是机械臂的重力矩阵，其具体形式分别如下：

由于现实工况中，考虑机械臂各个关节存在摩擦力f，则式（12）变为：

其中，摩擦力模型采用Stribeck模型，本小节考虑实际中关节正反向运动的摩擦力参数可能不同，因此摩擦力f形式如下：

其中，F_c、F_s、v_s、σ分别表示库伦摩擦力，静摩擦力，Stribeck速度，黏性摩擦系数.ζ是经验参数，（这里参考Lugre的静态模型，ζ=2）.此外，F_e表示外力，而正、负分别表示沿关节转动正方向和负方向.

本文假设Stribeck模型参数和机械臂各连杆的质量相关，参照文献^[32]中得到数据，可得Franka各关节摩擦力参数的信息，如下表：

2.2 基于Matlab的动力学仿真

本文以Franka七轴串联机械臂为对象，通过将Franka七轴串联机械臂模型导入到SolidWorks软件中，为模型施加材料属性（ρ=4650kg·m^-3）后测得其相关参数，见表3.

首先将表3的参数和式（15）代入到式（13）中，并通过Matlab对其进行计算得到机械臂的总动能T，接着将求得的T代入到式（17）和式（18），分别求得机械臂的质量矩阵和科氏矩阵，随后通过Matlab对式（14）进行计算得到机械臂的总势能U，最后将求得的U代入到式（19）中求得机械臂的重力矩阵，至此，基于Matlab的动力学方程建立初步完成.

2.3 动力学模型的验证

动力学模型的验证主要通过对比Matlab的动力学仿真输出的各关节轨迹图和基于Simulink和Adams的联合动力学仿真输出的各关节轨迹图.本小节所采用负载质量的大小为2kg，其体积为0.1m×1m×1m.

基于Matlab的动力学仿真本质上是一个正动力学问题，首先对驱动关节施加PD控制，即：

其中，

其中，q_e表示期望轨迹对应的关节旋转角度，q表示瞬时时刻的关节旋转角度.

令机械臂各关节的起始旋转角度均为0，而期望的关节旋转角度为：

接着采用三次插值法对各个关节进行轨迹规划，令t_d=5s，频率为1000Hz，即dt=0.001s，由此可得各关节角度的一组轨迹数据.

接着，将SolidWorks中已经附上材料属性的Franka七轴串联机械臂的模型导入到Adams软件中，并为各个关节设置约束条件和初始条件后导入到Simulink中添加Stribeck摩擦和PD控制（令T=5s，期望轨迹采用基于Matlab动力学仿真的期望轨迹，并利用ode45求解器进行求解），由此得到另一组轨迹数据.将Matlab得到数据与联合仿真得到的数据绘制成轨迹时程图进行比较，结果如图4所示.

从图4七个关节轨迹时程图对比来看，可以发现二者轨迹曲线定性吻合，因此验证了本文建立的动力学模型的有效性.

3 机械臂控制系统的稳定性分析

对于一个系统而言，无论其处于何种初始状态条件，如果其能量是耗散的，那么该系统在经过一个初始扰动后，就会不断损耗能量直至达到平衡点位置静止，这就是李雅普诺夫直接法的物理判据.该方法不需要通过求解复杂的系统非线性微分方程就可以得出系统是否稳定的结论，因此本小节采用李雅普诺夫直接法对动力学控制系统的稳定性进行分析^[33].

首先构造出李雅普诺夫函数，其形式如下：

因为质量矩阵D（q）和K_p都是正定的，因此该李雅普诺夫函数是正定的，当且仅当$\dot{\mathit{q}}$=0，e=0时，$V（\mathit{q}，\dot{\mathit{q}}）$=0.随后，对式（25）求一次导，因为定位控制有（${\dot{\mathit{q}}}_e$）≡0，所以可得：

结合式（16）、式（18）和式（26）可得：

令式（27）括号里的式子等于$-{\mathit{K}}_d\dot{\mathit{q}}$，即：

可得：

由于K_d是正定的，所以此时的$\dot{V}（\mathit{q}，\dot{\mathit{q}}）$≤0在整个空间中恒成立，当且仅当$\dot{\mathit{q}}$=0时，$\dot{V}（\mathit{q}，\dot{\mathit{q}}）$=0.所以式（28）控制律下的机械臂系统的定位控制能够实现全局渐近稳定.

随后对该控制律进行验证，采用P2（全局坐标为[-0.3  0.15  0.04]m）为期望位置，起始位置为P1（全局坐标为[-0.45  0  0.2]m），分别采用式（22）和式（28）的控制率，基于Matlab进行动力学仿真，并对其末端位置准确度分别进行测量，其末端位置准确度的计算公式如下：

其中，$\stackrel{-}{x}$、$\stackrel{-}{y}$、$\stackrel{-}{z}$表示多次对同一位姿重复响应后得到的位置平均值，而x_c、y_c、z_c表示指令位置.

经计算可得，式（22）控制下AP_p=8.24×10^-3 m，式（28）控制下AP_p=0.42×10^-7m.

随后对其末端的姿态准确度进行测量，其末端姿态准确度的计算公式如下：

其中，$\stackrel{-}{a}$、$\stackrel{-}{b}$、$\stackrel{-}{c}$表示多次对同一位姿重复响应后得到的姿态角的平均值，而a_c、b_c、c_c表示指令姿态角.

经计算可得，式（22）控制下AP_a=-0.0898、AP_b=-0.1781、AP_c=-6.5579×10^-4（单位弧度），式（28）控制下AP_a= 9.6593×10^-6、AP_b=-6.1366×10^-6、AP_c=4.6010×10^-5（单位弧度）.

对比位姿准确度，可以发现单独的PD控制并不能实现机械臂的定位控制，而添加重力补偿项和摩擦力补偿项的控制律能够很好地将定位控制的位姿误差维持在一个很小的范围内.

4 控制参数和负载质量对机械臂末端移动性能的影响规律

4.1 负载质量对机械臂末端移动性能的影响

本节采用定位控制，令起始位置为P1（全局坐标为[-0.4  0  0.2]m），期望位置为P2（全局坐标为[-0.3  0.15  0.04]m）.令末端负载质量由1kg每次增加0.5kg到5kg共9种不同的工况，进行动力学仿真，得到其中一末端位置准确度的时程图如下：

图5中越往外的曲线对应的负载质量越大，因此可以初步判断当末端质量越大时，机械臂末端误差收敛速率越慢.

为了验证该猜想，采用对数衰减率δ来体现其收敛速率v_AP，δ的定义如下：

其中A_i、A_i+1表示相邻两振幅.

由此可得末端负载质量与收敛速率的关系图.

观察图6，可发现末端误差收敛速率与质量之间的曲线呈现二次函数关系（v=am²+bm+c），因此本小节采用Matlab的Cftool工具对该曲线进行函数拟合，其拟合结果如下：

其中m表示末端负载质量，v_AP表示位置准确度误差收敛速率.

4.2 K_d参数和负载质量对机械臂末端移动性能的影响

本小节将考虑K_d参数、负载质量和末端收敛速度三者之间的关系.首先将K_d由2每次下降0.05到1.5，共11个不同控制参数取值（其中K_p≡15），最终得到误差收敛速率与末端质量的多项式函数各系数（a，b，c）与K_d之间的关系曲线图.

观察图7，可发现末端误差收敛速率各系数与K_d之间的曲线近似二次函数曲线（coefficient=$x{K}_d^2+y{K}_d+z$），其中，系数a与K_d成二次增加，系数b与K_d成二次减小，c与K_d成二次增加.因此本文采用Matlab的Cftool工具对该曲线进行函数拟合，其拟合结果如下：

从以上数据最终可以得到收敛速率v_AP与末端质量m与K_d之间的关系式为：

图8中，颜色越深表示数值越小，由此可得末端收敛速度与质量成负相关，与K_d值成正相关.

4.3 K_p参数和负载质量对机械臂末端移动性能的影响

本小节将考虑K_p参数、负载质量和末端收敛速度三者之间的关系.首先，将K_p由15每次增加0.5到20，共11个不同控制参数取值（其中K_d≡2），进行动力学仿真，最终得到误差收敛速率与末端质量的多项式函数各系数（a，b，c）与K_p之间的关系曲线图.

观察图9，可得末端误差收敛速率各系数与K_p之间的曲线近似二次函数曲线（coefficient=$x{K}_p^2+y{K}_p+z$），其中，系数a与K_p成二次减小，系数b与K_p成二次增加，c与K_p成二次减小.因此本文采用Matlab的Cftool工具对该曲线进行函数拟合，其拟合结果如下：

从以上数据最终可以得到收敛速率v_AP与末端质量m与K_p之间的关系式为：

图10中，颜色越深表示数值越小，由此可得末端收敛速度与质量成负相关，与K_p值成负相关.

5 结论

本文针对大负载引起的大惯性力不利于机械臂的提速增稳的问题，以Franka七轴串联机械臂为对象进行了运动学分析和动力学分析，随后对动力学控制系统的稳定性展开研究，采用李雅普诺夫直接法对定位控制的控制律进行设计并验证，结果表明，所设计的控制律能在满足系统需求的情况下实现全局渐近稳定.最后，基于定位控制对机械臂负载质量如何影响末端移动性能展开研究，其结果表明，负载质量增加，末端误差收敛速率减小（呈现二次函数关系）；当K_d值增加时，末端位置准确度收敛速率增大；当K_p值增加时，末端位置准确度收敛速率减小.

参考文献