文 | 猩猩说影

编辑 | 猩猩说影

机器人DH参数模型与静刚度模型

针对机器人进行运动学方向的建模分析是对其进行研究尤其是刚度辨识研究的必要条件，采用标准D-H法建立机器人的连杆坐标系，坐标系建立于结构中各杆后端的关节处。

且与该连杆固连，即前一坐标系i-1与对应的关节i对齐，同时规定该坐标系i-1的z轴即为机器人关节i的轴线，但其转向不完全跟该关节轴转向相同，限定坐标系x正向为相邻两z轴间的公共垂线方向，y轴方向则通过右手定则判断。

最终可得到的机器人各关节连杆坐标系，如图1所示，并把整个机器人加工系统的基坐标系标记为0，另外，机器人各连杆长度尺寸数据，如图2所示。

根据辨识的COMAUSMART5NJ165-3.0机器人数据，其D-H模型参数，如表1所示。

表1中，θ为对应关节的旋转角，α为相邻轴线间所成角度，a为z轴间沿公法向的距离长度，即对应连杆尺寸，d为相邻x轴间距离，即相邻杆之间的距离。

为使整体建模过程得到简化，降低具体计算难度，采用传统刚度模型完成对机器人加工系统刚度模型的建立，机器人关节多存在易弹性变形且刚度较低的传动部件，由于重载机器臂刚度大，在允许载荷范围内均可将各关节处的变形合量认定为机器人的末端法兰处变形。

基于外载荷作用下引起的变形统计数据，驱动、传动等部件引发变形于总变形中占绝大部分，因此首先进行变形假设：视机器人连杆为不可变形的刚性杆，而将关节视为扭簧，并以常量矩阵来表示关节刚度。

因此由假设可将机器人关节刚度矩阵表示为[Kθ]=diag[kθ1kθ2kθ3kθ4kθ5kθ6]的形式，式中：Ki—机器人第i个关节的刚度，且为常量也是待求量，机器人笛卡尔刚度矩阵K与其关节刚度矩阵Kθ的关系可表示为：

即为上文提到的传统刚度模型，式中：J—雅克比矩阵。

基于MATLAB的实验条件参数优化

针对刚度辨识试验设计加工了安装于机器人法兰处的专用加载测量装置，作为负载力施加位置，并进行变形测量，其中，装置留有激光跟踪仪靶球固定点，变形测量通过激光跟踪仪设备主体实现，因此结合装置形状尺寸以及试验现场的环境限制，需要在施力机构位置优化建模过程中注意以下几个问题：

负载力矢量不能与末端装置、机器人本体，现场其他物品和设备干涉，末端测量点需要易于测量，靶球反射不能偏离激光光轴，获得施力机构位置可以稳定放置，不能超出场地范围。

基于以上问题，分别在优化计算中进行了限制。

针对防止干涉的问题，从机器人末端法兰面开始分别将安装板、末端装置、激光跟踪仪靶球等各部分结构建立为圆柱模型，并表示负载力矢量，基于所建几何模型并通过函数求解力矢量与所有末端结构模型间是否存在交点，最终保留无交点结果。

另外限制施力装置放置范围，防止与机器人本体、现场其他物品和设备干涉，针对可测量性的问题，表示测量设备信号接收端点与机器人末端测量点间矢量，并求解各矢量与机器人第6轴轴线间的夹角，将所有夹角最大值限制于末端靶球可测量角度60°之内，从而保证可测量性。

因此，可将该位置优化问题视为有约束优化问题，并基于以上限制，在MATLAB程序中分别设置对应条件与取值范围，进行坐标与位姿角度计算，选用MATLAB中优化功能函数fmincon进行计算。

其中，计算所用的坐标数据均设置或转换为系统基坐标系下结果，考虑计算结果精度与时间合理性，对函数选项进行设置，算法选用序列二次规划法，并设置函数及变量的终止限度均为-1610，函数最大评价次数为20000。

在优化过程中，进行了多次不同初值选取，多个最终优化结果对比表明其结果与该因素无紧密联系，获取初始优化位置结果后，将施力结构摆放到对应位置以验证是否可用，从而保证了最终位置结果的可用性，多次迭代计算后，得到的最佳位置结果与机器人位姿组合，如表2与表3所示。

关节刚度辨识试验

试验设备包括COMAU机器人，FARO三维激光跟踪仪，负载施力机构，ATI六维力传感器及专用末端测量装置，施力机构通过钢丝绳将负载力施加到机器人末端。

激光跟踪仪通过设备与靶球测量机器人末端x、y、z三向位移变化，ATI力传感器连接法兰与专用测量装置并实时测量六自由度负载力信号与转矩信号，机器人型号为COMAUSMART5NJ165-3.0，其设计负载为165kg。

基于胡克定律分析机器人刚度计算中力与变形关系设计试验方案，其试验原理图如图3所示。

将测量装置安装于末端，调整机器人关节角度至预定位姿并将施力装置安放于优化后坐标，保证负载力的方向，向施力机构吊盘中添加负载，滑轮改变力的方向，实现将拉力施加到机器人末端。

试验时每次添加10kg重物负载，逐步添加至100kg，空载状态及每次添加负载后分别测量4个测量点对应变化的坐标5次，取其平均值作为当前位置坐标。

在测量靶球坐标的同时，记录力传感器示数并取平均值，改为基坐标系下实际负载分力，作为计算负载力的分力数据。

坐标值与分力值分别减去空载状态的值则可获得每组负载下的变形量与负载分力，基于机器人笛卡尔刚度模型中的力与变形规律，并根据已建立的机器人刚度模型，计算求出机器人的关节刚度，试验现场如图4所示。

试验结果数据基于机器人基坐标系，机器人笛卡尔刚度矩阵、负载力与变形量满足胡克定律F=KΔX，其中力矢量表示为：

变形矢量表示为：

式中：K-6×6机器人笛卡尔刚度矩阵，将式(1)代入胡克定律可得式(2)：

式中：[Kθ]-1—关节柔度矩阵[Kθ]C，即为关节刚度矩阵取逆表示为[Kθ]C=diag[kθ1-1kθ2-1kθ3-1kθ4-1kθ5-1kθ6-1]。

将[Kθ]C分离可得变形量ΔX与关节柔度矩阵[Kθ]C之间的关系：

式(3)中A矩阵仅与雅克比矩阵J与外力F有关，为便于计算，将[Kθ]C以6维向量的形式表示即为H=[kθ1-1kθ2-1kθ3-1kθ4-1kθ5-1kθ6-1]T，则A可表示为：

优化后的加载位置共有6个，且4个测量点位置均进行10组测量，并忽略转角和扭矩的影响，因此A矩阵最终维数为720×6，记作A0，变形矢量最终为720×1，记作ΔX0，由于矩阵与变形向量的规格问题，不能用正常矩阵求解手段，因此通过近似误差最小的近似解表示最终结果，即为式(5)：

通过矩阵伪逆计算，可获得误差最小近似[Kθ]C为

式中：A0l—A0的广义逆，因此机器人关节刚度矩阵[Kθ]=[Kθ]C-1，进行刚度计算，根据测量数据计算得到的刚度结果为(单位：N-mm/rad)

刚度计算过程中，由于获取数据不包括三方向的扭转分量，为确保结果的正确，广义力矢量的取值考虑了扭矩产生的影响，其数据进行了补偿。

在进行刚度辨识试验并求解获得关节刚度矩阵后，设计结果验证试验并记录数据，用以证明所求刚度矩阵可靠有效，根据所获刚度矩阵求解理论变形量，并计算所有实际变形量与理论变形量的差值。

选取4组可用位姿，在原有试验方案条件的基础上，分别添加50kg与100kg负载，同时4个测量点均记录变形量与力传感器示数，将验证试验获得的负载力利用所求刚度矩阵计算其对应变形量。

还选用了六组未优化位姿进行同样的刚度矩阵计算，并代入验证试验负载力获得对应理论变形量，在计算获得两组刚度对应理论变形后，分别求解其与验证试验测量得到的实际变形的偏差，并将差值统一取正可得两组偏差绝对值，如图5所示。

验证试验及数据对比结果证明通过该刚度辨识方法所得结果的准确性，表明提出方法切实可行，能够精确辨识的机器人系统的关节刚度矩阵，且相较任意位姿求解结果精确度有较大提升。

作者观点

通过综合优化施力装置的位置与机器人位姿的试验方式，而不是采用以往常用的简单挑选变换机器人位姿的手段，实现机器人刚度辨识，该方法相对以往试验方式，考虑因素更全面，其流程更为规范标准，为试验优化方向提供新的选择。

合理设置优化过程中的位姿角度取值范围，可有效避免位姿调整过程中的奇异位姿，也可保证了机器人灵巧性指标，因此，通过提出的刚度辨识试验方式，在保证结果准确的前提下也降低了风险，可获得更为合理的试验数据结果。

在刚度辨识试验中除位姿优化外，引入了负载施力装置安放位置优化，基于MATLAB软件为负载位置选择提供了规范化的程序，通过几何建模约束防止负载钢丝绳与其他试验设备及现场环境干涉。

在保证了可加载性的同时也考虑了变形可测量性，因此实现了试验方案及参数的标准化选择，其刚度辨识试验结果也更为接近实际结果，相对于常用的小负载辨识，在试验中选用了较大的负载，更加贴近实际加工工况，因此数据更加反映了真实结果。

提出的机器人刚度矩阵辨识方法，为机器人刚度分析提供了全新的研究手段，也为将来机器人加工研究提供了更为精确的刚度结果。