一种基于四阶矩估计的工业机器人可靠性计算方法与流程

本发明涉及了一种机器人可靠性计算方法，尤其是涉及了一种基于四阶矩估计的工业机器人可靠性计算方法，涉及工业生产、医药制造等领域中广泛使用的串联及并联型机器人。

背景技术：

随着工业4.0时代的到来，许多领域诸如大型复杂零部件的加工，装配和搬运，无夹具系统的焊接，喷涂以及雕刻等，对机器人的运动精度、稳定性及可靠性的要求也越来越高，但是由于机器人在设计、零部件加工、装配以及使用中难免磨损和受力变形等，使得机器人的本体存在着不可避免的误差，主要包括旋转关节处的关节间隙以及连杆尺寸的偏差，连杆尺寸的偏差使得机器人的末端运动精度降低，关节间隙的存在不仅影响运动精度，还会造成关节冲击，加剧振动甚至引起机器人的运动失效，因此对机器人在上述不确定性影响下的运动的可靠性研究具有重要意义。当前的工业生产中，用来对机器人的可靠性进行计算的方法主要包括了一阶可靠性分析方法、蒙特卡洛随机模拟方法，尤其是蒙特卡罗随机模拟方法更是广泛使用，但是由于蒙特卡洛模拟需要基于非常大的样本值才能较为准确的估计出系统的可靠性，在实际运用中，不仅成本高昂，而且效率也较低，一阶可靠性方法虽然可以减少模拟的样本数量，但是由于其计算精度有限，主要用在一些结构简单，非线性化程度较低的机械系统中，无法满足诸如航空航天、机器人运动控制系统等领域的要求。

技术实现要素：

为解决机器人可靠性计算的问题，本发明的目的是提出一种基于四阶矩估计的工业机器人可靠性计算方法，针对串联、并联以及任何旋转关节组成的机器人或机械系统，通过对关节间隙和连杆尺寸偏差的组成的随机变量进行统计学分析，确定其概率分布类型，通过机器人运动学获得机器人末端位置的误差组成的样本，建立了机器人系统的功能函数，并通过对末端位置误差进行统计分析，获得功能函数的概率密度函数表达式，从而快速准确的对机器人可靠性进行计算。

本发明采用的技术方案是采用以下步骤：

第一步：根据机器人的连杆臂结构建立包含基坐标系和关节局部坐标系的d-h连杆坐标系；

如附图2所示，首先将机器人的底部中心作为原点建立基座标系o_xyz，其中z轴沿机器人第一个旋转关节的轴向方向，x轴方向任意选取，y轴的指向根据右手法则由x轴和z轴确定。

如附图3所示，在第i个连杆臂旋转关节处建立一个局部坐标系oi-1xi-1yi-1zi-1，局部坐标系的原点oi-1为第i关节的旋转中心，zi-1轴沿第i关节的轴向方向，xi轴为当前关节轴向方向，yi轴的指向根据右手法则由xi轴和zi轴确定，则第i个连杆臂的局部坐标oixiyizi相对于第i-1个连杆臂的局部坐标之间的位置关系可通过如下变换得到：将局部坐标系oi-1xi-1yi-1zi-1绕zi-1轴旋转角度θi使得xi-1与xi轴平行，并沿zi-1轴平移di得到坐标系o′ix′iy′iz′i，将坐标系o′ix′iy′iz′i沿x′i平移ai并绕x′i轴旋转角度αi得到坐标系oixiyizi。

上述变换关系表示成一个4×4的齐次矩阵具体是指如下表达式：

其中，ai为第i个连杆臂长度，di为第i个连杆偏矩，αi为第i个连杆扭转角，θi为第i个连杆臂的关节转角，所述的机器人是串联或者并联的机械臂结构。

第二步：通过机器人连杆臂的正运动学获得机器人末端位置相对于基坐标系的齐次矩阵得到机器人末端的理论位置相对于基坐标系o_xyz的坐标值为p(x，y，z)，上述齐次矩阵采用如下表达式：

其中，n表示机器人的自由度数，r表示机器人末端在基座标系中的姿态矩阵，p(x，y，z)表示机器人末端位置坐标，п表示连乘符号；

第三步：将机器人第i个连杆臂处的连杆臂长度ai和连杆偏距di的尺寸偏差定义为满足正态分布的随机变量，连杆臂长度和连杆偏距的均值分别为μai和μdi，标准差分别为σai和σdi；同时将机器人第i个连杆臂处包含关节间隙的旋转关节的关节转角θi定义为满足均匀分布的随机变量，关节转角在均匀分布下的上边界为θi+δi，下边界为θi-δi；

第四步：利用matlab的随机数函数normrnd生成满足正态分布的随机数表示连杆臂长度和连杆偏距，利用随机数函数unifrnd生成满足均匀分布的随机数表示关节转角；

第五步：考虑机器人进行往复运动，在第j次往复运动中，由于连杆臂长度偏差、连杆偏距偏差和关节间隙的影响，机器人实际的末端位置会偏离理论位置，如附图3所示，偏离的距离定义为机器人末端运动误差ε，通过第四步生成的随机数值，利用第二步正运动学表达式的齐次矩阵计算机器人在第j次运动中实际末端位置相对于基座标系的坐标值为p′(x′，y′，z′)，再采用以下公式计算获得机器人末端运动误差εj：

第六步：设定机器人末端的定位精度为δr，以机器人的末端运动误差ε小于定位精度δr的概率作为机器人运动可靠性，将实际机器人末端运动误差ε视为随机变量，机器人每进行一次往复运动，则存在一个末端运动误差值，重复第四步～第五步进行n次，共获得n个机器人末端运动误差值，将所有机器人末端运动误差值组成一个样本数据；

第七步：对第六步中获得的样本数据进行分析，计算机器人末端运动误差ε的均值με，计算公式为：

其中，εj表示第j次运动中的机器人末端运动误差；

第八步：对第六步中获得的样本数据进行分析，计算机器人末端运动误差ε的前四阶中心矩计算公式为：

其中，k表示阶数，取值为k＝0，1，2，3，4；

同时以二阶中心矩作为机器人末端运动误差的方差，即其中σε表示机器人末端运动误差的标准差；

第九步：根据机器人运动可靠性的定义，建立机器人系统的功能函数g(ε)，表达式为：

z＝g(ε)＝δr-ε

其中，z表示功能函数的取值，δr表示机器人末端的定位精度，ε表示机器人末端运动误差；

第十步：根据机器人末端运动误差ε的均值με及其前四阶中心矩和通过第九步建立的机器人系统的功能函数g(ε)，计算机器人系统的功能函数g(ε)的均值μz，计算的公式为：

μz＝δr-με

再计算机器人系统功能函数的前四阶中心矩计算公式为：

其中，k表示具体的阶数，取值为k＝0，1，2，3，4；

机器人系统的功能函数的零阶中心矩一阶中心矩二阶中心矩三阶中心矩四阶中心矩

同时以二阶中心矩作为机器人系统功能函数的方差，即其中σz表示机器人系统功能函数的标准差；

第十一步：将功能函数值z变换为均值为0、方差为1的标准随机变量y，变换公式为：

计算标准随机变量y的前四阶中心矩计算公式为：

其中，k表示具体的阶数，取值为k＝0，1，2，3，4；

第十二步：定义随机变量y的信息熵为h[f(y)]，信息熵具体表达式为：

h[f(y)]＝-∫yf(y)ln[f(y)]dy

其中，f(y)表示随机变量y的概率密度函数；

将标准随机变量y的前四阶中心矩作为约束条件，根据最大熵原理构造拉格朗日函数，拉格朗日函数表达式为：

其中，λk为第k阶对应的拉格朗日系数，λk＝(λ0，λ1，…，λ4)；

利用拉格朗日函数l对概率密度函数f(y)求偏导数，使得在极值点处偏导数为0，求导过程具体表达式为：

第十三步：根据上一步求导过程，采用以下公式计算获得随机变量y的概率密度函数f(y)，表达式具体为：

第十四步：由于机器人的可靠性是指机器人在连杆尺寸偏差、连杆偏距偏差以及关节间隙的影响下，机器人末端的运动误差ε在给定的定位精度δr范围之内的概率，即在某一次运动过程中，机器人的末端运动误差ε<δr的概率，所以经功能函数后，机器人运动的可靠性进一步转换成：

其中，pr{z≤0}表示满足z≤0条件的概率，pr{z≤0}等价于

根据随机变量y的概率密度函数f(y)，计算出机器人的可靠性计算过程表达式具体为：

其中，σz表示机器人系统功能函数的方差，μz表示机器人系统的功能函数g(ε)的均值。

所述第二步中，包括关节间隙和连杆尺寸偏差，连杆尺寸主要包括机器人dh参数中的连杆长度和连杆偏距，因此连杆尺寸偏差包括了连杆长度偏差和连杆偏距偏差。由于连杆臂在加工制造、装配过程中不可避免存在误差，主要表现为连杆尺寸与理论值之间存在差别，因此所述误差参数包括关节间隙和连杆尺寸偏差，连杆尺寸主要包括机器人dh参数中的连杆长度和连杆偏距，因此连杆尺寸偏差包括了连杆长度偏差和连杆偏距偏差，机器人的连杆长度偏差及连杆偏距偏差主要是因为加工、装配以及受力变形等因素引起，因此连杆长度偏差是满足正态分布的随机变量。

所述第三步中：如附图4所示，机器人相邻的连杆臂之间是通过转轴和轴承组成的旋转副连接而成，因此关节间隙表现为转轴和轴承的旋转中心不重合，在机器人运行过程中，转轴在轴承的内部随机窜动，造成关机冲击和振动，使得关节的实际转角和理论转角之间存在误差，因此以转轴和轴承的旋转中心之间的距离作为关节间隙，转轴的旋转中心落在轴承内圈构成的一个圆形区域内，落点位置是随机的并且是等概率的，即关节间隙的误差是符合均匀分布的随机变量。

所述第四步中：本发明中利用matlab随机数函数normrnd生成了包含偏差并且符合正态分布的连杆长度参数和连杆偏距参数，利用随机函数unifrnd生成了包含关节间隙并且符合均匀分布的关节转角参数，本发明方法中，机器人的所有连杆臂均存在连杆尺寸误差，所有旋转关节均存在关节间隙。

所述第五步中，不考虑关节间隙和连杆尺寸偏差的影响下，机器人末端通过正学运动学方程计算的实际位置和理论位置完全重合，但由于机器人的关节处存在关节间隙以及连杆尺寸存在一定误差，由此使得机器人实际末端位置和机器人理论末端位置并不重合，存在一定的误差，因此以实际位置偏离理论位置的距离作为机器人的运动误差。

所述第六步中，机器人的关节间隙和连杆尺寸偏差分别是满足均匀分布和正态分布的随机变量，因此在运动中，机器人的各个关节转角和连杆尺寸参数是存在误差的，通过正运动学方程获得的机器人末端实际位置和理论位置之间存在误差，由于关节间隙和连杆尺寸偏差的随机性，机器人的运动误差也具有随机性，所以将其视为随机变量。

所述第七步中，由于机器人末端的运动误差具有随机特性，因此使得功能函数值也具有一定的随机特性，所以将功能函数值视为随机变量。因此通过统计学方法，可以对由机器人末端的运动误差组成的样本数据进行分析，分别获得表征机器人末端的运动误差统计特性的一阶原点矩(均值)με。

所述第八步中，计算出机器人末端运动误差的零阶中心矩一阶中心矩二阶中心矩(方差)三阶中心矩以及四阶中心矩

所述第十步中，然后通过功能函数表达式，利用机器人末端运动误差的均值进一步获得功能函数值的一阶原点矩(均值)μz，零阶中心矩一阶中心矩二阶中心矩(方差)以及三阶中心矩和四阶中心矩

所述第十一步中，为了简化计算过程，将机器人功能函数z转换成标准随机变量y，并计算了标准随机变量y的均值μy，零阶中心矩一阶中心矩二阶中心矩(方差)以及三阶中心矩和四阶中心矩

所述第十二步中，机器人的运动可靠性取决于标准随机变量y的概率分布，因此为了计算可靠性，需要获得标准随机变量y的概率密度函数模型表达式，标准随机变量y的均值以及前四阶中心矩已经可以得到，利用概率论的最大熵原理，构造拉格朗日函数，后续进一步将标准随机变量y的概率密度函数模型的求解问题转换成满足一定约束条件的优化问题。

所述第十三步中，通过求解满足约束条件下的优化问题获得拉格朗日函数中的未知参数，从而得到完整的具有参数的标准随机变量y的概率密度函数模型表达式，通过对概率密度函数模型的积分，实现机器人的运动可靠性计算。

本发明利用了关节间隙及连杆长度偏差及连杆偏距偏差的概率分布模型，通过机器人的正运动学建立了机器人系统的功能线性函数，避免了对机器人的末端位置运动函数的泰勒级数展开，从而消除了由于将非线性函数线性化带来的误差，同时，因为关节间隙、连杆长度偏差及连杆偏距偏差的随机性使得机器人末端位置误差也具有一定的随机性，本发明将机器人末端运动误差视为随机变量，利用统计分析求取了其均值和前四阶中心矩，充分利用了随机变量的统计特征，最大程度上避免了对随机变量包含的统计信息的丢失，利用末端位置误差的统计特征进一步计算出了功能函数的相关统计特征，并基于最大熵原理求取了功能函数的概率密度函数表达式，最后利用积分，获得机器人的可靠性。

本发明的优势和有益效果在于如下几点：

1、针对含关节间隙、连杆长度偏差及连杆偏距偏差的机器人，本发明解决了机器人可靠性分析问题，作为机器人设计、加工、装配的前提，是保证机器人在实际运用之中能够精确、稳定、可靠运行的技术基础。

本发明在计算过程之中，首先基于关节间隙、连杆长度偏差以及连杆偏距偏差的概率模型分析了机器人末端运动误差的统计特征，将末端运动误差视为随机变量，根据机器人定位精度及可靠性定义，建立了机器人系统的线性功能函数，避免了直接对机器人的正运动学方程进行泰勒级数展开，避免了将非线性函数进行线性化所引起的误差。

2、本发明方法，在可靠性分析中相较于目前广泛使用的蒙特卡洛模拟法，在保证结果准确性的前提下，能够大幅度降低对样本规模的需求，仅需极少的样本数量便获得了较高的计算精度。

3、本发明方法，在可靠性分析中，相较于传统的一阶可靠性分析法，能够利用随机变量更多的统计信息。

传统的一阶可靠性分析法仅仅利用了变量的一阶原点矩(均值)和二阶中心矩(方差)的统计信息，在变量不严格满足正态分布情况下会存在较大误差。

本发明方法不仅利用了随机变量的均值和方差，还进一步利用了随机变量的三阶(反映随机变量概率密度函数的偏度信息)及四阶中心矩(反映随机变量概率函数的峰度信息)的统计信息，虽然两者所需的样本数量均较少，但是本发明方法能够获得更加准确的计算结果。

附图说明

图1是本发明方法的流程框图；

图2是本发明中适用的旋转关节型机器人示意图；

图3是本发明中机器人的连杆坐标系示意图；

图4是本发明中机器人末端运动误差示意图；

图5是本发明中机器人旋转关节间隙示意图。

具体实施方式

下面结合附图及具体实施例对本发明做进一步详细说明

本发明的基本原理主要在于：机器人的关节间隙是满足均匀分布的随机变量，连杆长度偏差及连杆偏距偏差是满足正态分布的随机变量，关节间隙会造成机器人的关节转角存在一个在一定范围内均匀分布的随机误差，连杆尺寸偏差会使得机器人的物理参数发生变化，关节间隙、连杆长度偏差及连杆偏距偏差通过机器人的运动学会造成机器人末端实际的位置偏离其预期位置，实际位置和理想位置之间的距离是机器人末端位置的运动误差，该误差也属于随机变量，根据机器人的末端定位精度，建立了机器人系统的功能函数，通过功能函数建立了机器人可靠性的计算公式，功能函数的取值也属于随机变量，利用由机器人末端运动误差组成的样本，进行统计学分析，得到表示运动误差的随机变量的均值以及前四阶中心矩，利用其均值和前四阶中心矩通过功能函数进一步求取功能函数的均值和前四阶中心矩，将功能函数的前四阶原点矩作为约束条件，基于概率论的最大熵原理，将该问题进一步转换成了满足约束条件下的优化求解问题，建立拉格朗日方程，求解其拉格朗日系数，从而得到功能函数的概率密度函数，通过积分，计算出机器人的可靠性。

本发明的实施例：

附图1中的框图为本发明提出的新型机器人可靠性计算方法的具体流程图，如附图所示，如附图2所示，本发明以一台3自由度串联型工业机器人为分析对象，包括连杆臂1和旋转关节2，针对含关节间隙及、连杆长度偏差及连杆偏距偏差的机器人可靠性计算，结合实施例进行说明，其中涉及的角度均以度(°)表示，涉及的尺寸长度以及距离的单位均为毫米(mm)，本发明的具体流程主要包括以下步骤：

第一步：

建立机器人的d-h连杆坐标系，机器人相应的物理参数包括如下：连杆长度a2＝400(mm)，a3＝600(mm)，连杆偏距d1＝500(mm)，在目标点处，机器人的3个旋转关节的理论转角分别为：θ1＝100°，θ2＝-60°，θ3＝50°，则机器人末端的理论位置坐标为p(x，y，z)，其中x＝-137.59(mm)，y＝780.36(mm)，z＝-206.43(mm)。

第二步：

将连杆臂自身和连杆臂之间的关节间隙、连杆长度偏差及连杆偏距偏差定义为满足一定分布类型的随机变量，在误差参数的范围内随机产生满足分布类型的随机数值作为误差参数，如图4所示；

关节间隙具体举例说明：具体实施的两个关节连接处如图5所示，一个关节的连杆3端部与转轴1固定连接，另一个关节的连杆3端部与轴承2固定连接，转轴1套装在轴承2中，由于转轴1外圈和轴承2内圈存在间隙，因此造成了转轴和轴承的旋转中心不重合，并以转轴和轴承的旋转中心之间的距离作为中心偏距4，使得关节在旋转过程中存在一个关节间隙δ。

连杆尺寸偏差属于正态分布，各自的均值和标准差分别为：

μa2＝400(mm)，σa2＝0.35(mm)，μa3＝600(mm)，

σa3＝0.30(mm)，μd1＝500(mm)，σd1＝0.35(mm)

关节间隙δ满足均匀分布，范围为[-0.1°，0.1°]，则机器人的实际关节转角为：

θ′i＝θi+δ(i＝1，2，3)

在连杆尺寸偏差及关节间隙的影响下，机器人在第j次往复运动中，末端实际位置坐标为pj′(xj′，yj′，zj′)，则机器人末端位置误差为：

第三步：

重复步骤二500次，获得由500个由末端运动误差组成的随机样本，进行统计学分析，分别可以求得其均值με＝1.284(mm)，零阶中心矩一阶中心矩二阶中心矩(方差)三阶中心矩四阶中心矩

第四步：

给定机器人末端的定位精度为δr＝2(mm)，由此建立机器人的功能函数z＝g(ε)＝δr-ε，利用随机变量ε的均值及前四阶中心矩计算出功能函数z的均值及前四阶中心矩分别为：μz＝0.716(mm)，根据机器人可靠性的定义：在关节间隙及连杆尺寸偏差影响下，机器人的末端位置误差小于定位精度的概率，所以机器人的可靠性可以进一步表达成：

第五步：

将功能函数z转换成均值为0，方差为1的标准随机变量y，转换关系可以表示成：

根据z的前四阶中心矩可以求得y的前四阶中心矩，对应关系为：

将z的相关数据代入可以分别求得y的前四阶中心矩为：因为随机变量y为标准随机变量，即均值为0，方差为1，所以其前四阶中心矩等于前四阶原点矩。

第六步：

基于概率论中最大熵原理，将随机变量y的前四阶原点矩作为约束条件，并定义

随机变量y的信息熵为：

h[f(y)]＝-∫yf(y)ln[f(y)]dy

其中，f(y)为随机变量y的概率密度函数，

建立拉格朗日方程为：

根据优化求解过程得到随机变量y的概率密度函数估计：

其中，λ＝(λ0，λ1，…，λ4)为拉格朗日系数。

第七步：

根据第五步，随机变量y的前四阶中心矩已经求得，利用第六步建立的随机变量y的概率密度函数表达式，构造包含拉格朗日系数的方程组：

利用该方程组可以求解得到拉格朗日系数λ＝(λ0，λ1，…，λ4)，进而确定随机变量y的概率密度函数表达式fy(y)，将前述步骤中的数据代入，求解得到拉格朗日系数λ＝(λ0，λ1，…，λ4)分别为：λ0＝0.0437，λ1＝0.0617，λ2＝0.2935，λ3＝-0.1538，λ4＝1.0131，从而可以得到随机变量y的概率密度函数f(y)的表达式。

第八步：

根据第四步、第五步，机器人的可靠性可以转换成

根据第七步所求的随机变量y的概率密度函数，机器人可靠性可以通过下述方法求得：

在考虑关节间隙及连杆尺寸偏差影响下，给定末端定位精度，机器人的末端位置误差小于定位精度的概率为，即机器人的可靠性为

实施例以传统较权威的蒙特卡洛模拟法进行同样实施所获得的可靠性为0.91236。

由此可见，本发明方法能准确地计算获得工业机器人的可靠性，

并且具体实施中，传统的蒙特卡洛模拟法利用了10万个样本值，计算时间为1182.63s本发明利用的样本值为500个，计算时间为12.64s相较于传统的蒙特卡洛模拟法有大幅度的减少，效率上可以提高近100倍，并且计算精度高。

由此可见本发明不仅具有效率高、计算精度好的等优点，而且大大减少了计算量，减少了计算时间，而在机器人可靠性分析中具有显著技术效果。