基于Copula函数的工业机器人电气驱动器的可靠性分析方法与流程

本发明属于工业机器人主要部件的可靠性分析技术领域，特别涉及一种基于copula函数的工业机器人电气驱动器的可靠性分析方法。

背景技术：

工业机器人驱动器是机器人完成工作的动力装置，作用相当于人的肌肉，是决定工业机人作业性能的核心组成部件。其中电气驱动器具有惯量低、速度与位置精度高、速度调节范围大的优点，并且电气驱动器不需进行能量转换，使用控制方便、灵活，在工业机器人中被广泛使用。

随着制造业的发展，对于工业机器人运动速度、定位精度、承载能力等性能的要求不断提高，使得对电气驱动器的性能参数要求也不断提高，同时导致了电气驱动器的工作环境与工作载荷较以前更为复杂，如由于环境气温和运行发热导致的高温与温度循环，由于气候导致的高湿度，由于电网能量输入不稳定和负载不稳定导致的电压过高与大幅波动，从而使其失效问题愈来愈突出，一旦电气驱动器失效，必然导致整个工业机器人系统的故障，造成不可估计的经济财产损失，甚至威胁人身安全。因此，研究工业机器人电气驱动器的可靠性，具有重要的实际意义。

目前，针对电气驱动器的可靠性分析研究，往往只关注某单一的失效模式，或者仅仅把多种失效模式看作是独立的串并联系统，并不能真实而有效地反应电气驱动器各个模块组成部分之间的结构关系，对整个驱动器系统的可靠性评估，往往会得到过低或过高的结果。

技术实现要素：

本发明的目的在于克服现有技术的不足，提供一种使用copula函数的方法描述各个主要失效模式之间的相关关系，建立基于copula函数的工业机器人电气驱动器的可靠性分析方法，克服了因未考虑电气驱动器各主要失效模式之间的相关关系，仅将它们视为串联模型，而导致在可靠性分析中低估其可靠性的问题。同时在保证精度的前提下，提升了对于电气驱动器使用传统的montecarlo仿真方法进行可靠性分析的计算效率。

本发明的目的是通过以下技术方案来实现的：基于copula函数的工业机器人电气驱动器的可靠性分析方法，包括以下步骤：

s1、根据电气驱动器试验和维修的统计数据，将电气驱动器的各类失效模式进行汇总和分析，使用故障模式、影响及危害性分析方法确定出驱动器的主要失效模式，并找出与各主要失效模式所对应的元器件；

s2、建立与驱动器各主要失效模式相对应的极限状态方程；

s3、使用montecarlo法产生montecarlo仿真样本，带入各失效模式的极限状态方程中求得响应值组g1，g2，…，gm，gi＝gi(x)，其中gi(x)为第i个极限状态方程，x为montecarlo仿真样本，并计算各失效模式发生的概率；

s4、根据步骤s3中产生的montecarlo仿真样本，定性分析各失效模式之间的相关关系，确定备选的copula函数类型；

s5、根据步骤s3中产生的montecarlo仿真样本，从步骤s4给出的备选copula函数类型中选择最佳copula函数并确定函数的参数值；

s6、根据步骤s5所得出的copula函数，计算考虑多失效模式相关的工业机器人电气驱动器失效概率与可靠度。

进一步地，所述步骤s1中通过使用故障模式、影响及危害性分析方法确定出的电气驱动器主要失效模式包括：电源与功率驱动模块中，电容的电参数(容量和等效串联电阻)恶化失效；功率驱动模块中，igbt模块铝键合线脱落断裂失效；印制电路板中，镀通孔疲劳断裂失效。

进一步地，所述步骤s2中采用基于故障物理方法对电气驱动器的失效模式建立极限状态方程。具体包括以下子步骤：

s21、使用故障物理方法分析各失效模式的失效产生原因；

s22、根据加工工艺、环境因素和工作载荷情况，分析失效原因中的不确定因素；

s23、根据步骤s21与s22所知的失效模式的特点，确定其所需要使用的故障物理模型；

s24、根据步骤s22所知的不确定因素，量化步骤s23所得的故障物理模型中参数，确定随机变量的分布类型与分布参数，建立失效模式的极限状态方程，记为gi(xi)＝0，其中i为失效模式的序号，xi为关于主要失效模式i的随机变量向量。

进一步地，所述步骤s4中定性分析各失效模式之间的相关关系包括以下子步骤：

s41、分别求取步骤s3中产生的每组响应值g1，g2，…，gm的经验分布函数，记为其中，其中n为montecarlo仿真样本的个数，x为经验分布函数的自变量，gij为第j个montecarlo仿真样本带入第i个极限状态方程所求得的响应值，1gij≤x为事件gij≤x的指示函数，当gij≤x时1gij≤x＝1，否则1gij≤x＝0；

s42、将步骤s3中每组样本点的响应值组g1，g2，…，gm带入各自的经验分布函数中，得到新的响应值组，记为y1，y2，…，ym；

s43、两两绘制出步骤s42所得新的响应值组y1，y2，…，ym的散点图，观察分析每两种失效模式之间的相关关系，选出备选的copula函数类型：当散点图图像呈现椭圆分布的形状时，采用椭圆族copula函数(正态copula、t-copula等)；当散点图图像呈现上下尾相关的形状时，采用阿基米德族copula函数(gumbelcopula、claytoncopula和frankcopula等)；若散点图中点均匀分散地分布在图中的每处，则说明两种失效模式相关关系很小，不需采用copula函数进行分析，此时视为相互独立即可。

进一步地，所述步骤s5中使用kendall秩相关系数确定copula函数的参数值，并选择平方欧氏距离最小的copula函数作为最优函数。步骤s5包括以下子步骤：

s51、两两求取步骤s42所得的y1，y2，…，ym的kendall秩相关系数；

s52、通过步骤s51中所得的kendall秩相关系数，使用解析法分别求出每两种失效模式组的备选copula函数的参数值；

s53、对每组失效模式分别求出每一种备选copula函数的平方欧氏距离，选择平方欧氏距离最小的copula函数为最优copula函数。

本发明的有益效果是：

1、本发明确定了工业机器人在复杂工况下，其使用的电气驱动器的主要失效模式、薄弱环节以及失效原因，并针对主要失效模式对应的元器件建立了极限状态方程，比传统仅依靠统计数据进行电气驱动器可靠性分析更为精确。

2、本发明使用copula函数的方法描述了各个主要失效模式之间的相关关系，建立了基于copula函数的电气驱动器的可靠性分析方法，克服了因未考虑电气驱动器主要失效模式之间的相关关系，仅将它们视为串联模型，而导致在可靠性分析中低估其可靠性的问题。同时在保证精度的前提下，提升了对于电气驱动器使用传统的montecarlo仿真方法进行可靠性分析的计算效率。

附图说明

图1为本发明的可靠性分析方法的流程图；

图2为本发明的电气驱动器结构图；

图3为igbt模块铝键合线脱落断裂失效结构示意图；

图4为igbt模块铝键合线脱落断裂失效实际情况图；

图5为本发明实施例中，电容的电参数(容量与等效串联电阻)恶化失效的相关关系散点图；

图6为本发明实施例中，电容的电参数(容量)恶化与igbt模块铝键合线脱落断裂失效的相关关系散点图；

图7为本发明实施例中，电容的电参数(容量)恶化与镀通孔疲劳断裂失效的相关关系散点图；

图8为本发明实施例中，电容的电参数(等效串联电阻)恶化与igbt模块铝键合线脱落断裂失效的相关关系散点图；

图9为本发明实施例中，电容的电参数(等效串联电阻)恶化与镀通孔疲劳断裂失效的相关关系散点图；

图10为本发明实施例中，igbt模块铝键合线脱落断裂与镀通孔疲劳断裂失效的相关关系散点图。

具体实施方式

copula函数模型是基于sklar定理提出的一种新型的多元联合分布建模工具，它能捕获联合分布中各变量之间完整的相关关系，摆脱了边际分布对联合分布整体的影响，被认为是研究变量相关性方法中最合适的选择之一。目前copula函数模型无论是在金融、经济等社会科学领域，还是在水文、地质等自然科学领域都得到了广泛的应用，同时其也已经在可靠性研究方面受到了特别的重视。所以使用copula函数模型来处理本发明中的对象和相关问题，更具有意义。

下面结合附图进一步说明本发明的技术方案。

如图1所示，一种基于copula函数的工业机器人电气驱动器的可靠性分析方法，包括以下步骤：

s1、根据电气驱动器试验和维修的统计数据，将电气驱动器的各类失效模式进行汇总和分析，使用故障模式、影响及危害性分析方法确定出驱动器的主要失效模式，并找出与各主要失效模式所对应的元器件；

电气驱动器，一般可以分为：电源、功率驱动、保护、控制、通信、交互、检测等模块。关系结构图如图2所示，其中电源、功率驱动、控制、检测模块较为重要：电源模块为电气驱动器中其他各硬件部分提供能量；功率驱动模块对伺服电机提供能量，并控制电机转速转向；控制模块通过通信模块与控制器相连接，接收控制器的命令，通过检测模块获得电机实时状态，通过功率驱动模块对电机进行控制。

在本实施例中，根据电气驱动器试验和维修的统计数据，使用故障模式、影响及危害性分析方法选择风险优先数大于200的故障模式，可以得到电气驱动器主要失效模式包括：电源与功率驱动模块中，电容的电参数(容量与等效串联电阻)恶化；功率驱动模块中，igbt模块铝键合线脱落断裂；印制电路板中，镀通孔疲劳断裂。

s2、根据步骤s1中所确定的驱动器各主要失效模式，建立与驱动器各主要失效模式相对应的极限状态方程；采用基于故障物理方法对电气驱动器的失效模式建立极限状态方程，具体包括以下子步骤：

s21、使用故障物理方法分析各失效模式的产生原因；

s22、根据加工工艺、环境因素和工作载荷等情况，分析失效原因中的不确定因素；

s23、根据步骤s21与s22所知的失效模式的特点，确定其所需要使用的故障物理模型；

s24、根据步骤s22所知的不确定因素，量化步骤s23所得的故障物理模型中参数，确定随机变量的分布类型与分布参数，建立失效模式的极限状态方程，记为gi(xi)＝0，其中i为失效模式的序号，xi为关于主要失效模式i的随机变量向量。

本实施例建立的主要失效模式的极限状态方程分别为：

(1)铝电解电容的电参数(容量)恶化

在本实施例中，铝电解电容的容量参数恶化的表现形式为电容量的下降到其所能被接受的阈值(标称电容量的80％)。该失效模式的失效原因为在工作过程中，铝电解电容的阴阳两极会发生化学反应(阳极：2al+3o^2-＝al2o3↓+6e^-，阴极：2h⁺+2e^-＝h2↑)，使得阳极的金属氧化膜不断被修补增厚，导致电容量c不断下降；失效的影响因素主要有环境的温度、湿度，工作的电压，电容的生产工艺、尺寸等。

根据以往的恒定应力加速寿命实验所得的数据，当退化速率υca＝ddc/dt恒定时，电容量c(t)与时间t呈线性下降关系，即

c(t)＝c(0)×(1-υcat)(1)

而在实际工作中，退化速率υ＝dd/dt随温度、湿度、工作电压等因素变化，所以采用广义艾琳(eyring)加速模型进行退化速率的描述，可以得到实际退化速率υc与恒定应力加速寿命实验所求得的退化速率υca的比值γc

其中，为退化反应激活能与玻尔兹曼常数的比值，t与ta分别为实际工作电容核心温度与额定工作温度，v与va分别为实际工作电压与额定工作电压，rh与rha分别为实际工作环境相对湿度与额定环境相对湿度，α与β分别为电压系数与湿度系数。

而电容在大批生产过程中，又因为材料、生产工艺等方面存在不一致性，电容出厂时的电容量c(0)会和标称电容量c0之间存在一定的误差，称为容量误差δc，那么电容量c(t)可以如下求得：

所以在t时刻的极限状态方程为：

(2)铝电解电容的电参数(等效串联电阻)恶化

铝电解电容在工作过程中，其电解质溶液中的水分会不断蒸发，使得其不断变稠。而电解液的电阻率会因电解液的粘稠度变大而增加，导致电容的等效串联电阻esr不断上升，最终超过所规定的阈值(额定值esr0的2.5倍)而使电容失效。

根据以往的恒定应力加速寿命实验所得的数据，当退化速率υesra＝dde/dt恒定时，等效串联电阻esr(t)与时间t呈指数上升关系，即

esr(t)＝esr(0)×exp(υesrat)(5)

同样的，在实际工作中的退化速率也可以使用艾琳模型进行描述，实际退化速率υesr和恒加实验求得的速率υesra比值为γesr；而又因为电解液(如硼酸——乙二醇)的粘稠度会在低温条件下升高，有所以等效串联电阻esr(t)为：

那么在t时刻的极限状态方程为：

(3)igbt模块铝键合线脱落断裂

当igbt模块长期工作在热循环冲击下，由于模块内材料的线性热膨胀系数不同，会产生热应力导致材料疲劳和老化，最终导致模块失效。在本实施例中，因为在igbt模块中芯片与焊料层的线性热膨胀系数之差为25×10^-6/k、铝键合线与芯片层的线性热膨胀系数之差为19×10^-6/k，较之其他层大，所以一般情况下失效多发生于此。igbt模块铝键合线脱落断裂失效结构示意图如图3所示，首先焊料层因为热应力导致的开裂会使得芯片工作产生的热量无法向下传递，从而热量聚集在铝键合线处，产生较大的热应力，而热应力又导致较大的剪切应力，从而使得键合线产生裂纹，最终脱落断裂失效；对于铝键合线脱落失效实际情况，如图4所示。

其中影响失效的因素有工作时温度的循环幅度，生产的工艺、尺寸等。而一般使用的故障物理模型为基于应变的manson-coffin模型，其通用公式为：

其中，nf为键合点发生破坏时其所经历的热循环次数，即疲劳失效循环次数；δεp为非弹性剪切应变变化幅度；εf为疲劳延性系数，取值与键合的材料组成成分有关，在本实施例中取εf＝0.257；c为疲劳延性指数，在本实施例中取c＝-0.415。而根据wernerengelmaier的研究发现剪切应变变化幅度δεp与具体电子部件的形式、几何尺寸、材料属性及其温度载荷等有关。所以最终极限状态方程为：

其中，ld为焊点的长度，h为焊点的高度，δt为周期性温度加载范围，δα为键合线与焊接层线性热膨胀系数的差值。

(4)镀通孔疲劳断裂

镀通孔的失效，主要是因为镀层材料与基板材料之间的线性热膨胀系数不匹配造成的，主要表现在印制电路板的厚度方向和最外层的焊盘上，最外层的焊盘的线性热膨胀系数通常为镀层线性热膨胀系数的3到4倍。当印制电路板在整个寿命周期内温度环境条件十分复杂的情况下(如制造和使用过程中的热熔、焊接以及电路工作发热导致的温度波动)，线性热膨胀系数的不匹配会使镀通孔中产生热应力导致镀层的疲劳损伤，并最终失效。

其中影响失效的因素主要是镀通孔加工的形状尺寸与位置。而镀通孔的疲劳断裂失效，本质上是镀层铜箔的疲劳断裂失效，所以使用一系列铜箔实验得出的实验结果模型为：

其中，nf为失效前平均循环周期数，df为pth镀层材料的断裂应变，su为pth镀层材料的断裂强度，ecu为金属(铜)镀层的弹性模量。

由应力模型可以分析得到，镀通孔镀层所受轴向正应力在中心处最大，即最大应变发生在中心处，所以用此处的应变δε带入上述模型，即可得到镀通孔疲劳断裂失效的极限状态方程为：

其中，为树脂基板厚度的一半；αcu、αe分别为金属(铜)镀层和树脂基板材料的线性热膨胀系数；δt为温度循环幅值；r为基板有效作用半径(即镀通孔中心到剪应力τps＝0处的距离)；为镀通孔的半径；t为镀通孔镀层的厚度；ge为树脂基板材料的剪切模量；k为应力校准系数，一般取k＝0.25。

s3、使用montecarlo法产生montecarlo仿真样本，带入各自的极限状态方程中求得响应值组g1，g2，…，gm，gi＝gi(x)，其中gi(x)为第i个极限状态方程，x为montecarlo仿真样本，并计算各失效模式发生的概率；在本实施例中，使用matlab软件进行montecarlo仿真，生成n＝1×10⁶个随机的多维样本点，并带入步骤s2所得到的极限状态方程中，求出响应值组g1，g2，g3，g4：g1代表电容的电参数(容量)恶化失效的响应值组，g2代表电容的电参数(等效串联电阻)恶化失效的响应值组，g3代表igbt模块铝键合线脱落断裂失效的响应值组，g4代表镀通孔疲劳断裂失效的响应值组。分别对响应值组g1，g2，g3，g4中小于0的样本进行计数，除以样本点的个数n＝1×10⁶，可得每种主要失效模式发生的概率：电容的电参数(容量)恶化发生的概率为pf1＝0.0192，电容的电参数(等效串联电阻)恶化发生的概率为pf2＝0.0113，igbt模块铝键合线脱落断裂发生的概率为pf3＝0.0034，镀通孔疲劳断裂发生的概率为pf4＝0.0193。

s4、根据步骤s3中产生的montecarlo仿真样本，定性分析各失效模式之间的相关关系，确定备选的copula函数类型；所述定性分析各失效模式之间的相关关系使用matlab软件进行分析，包括以下子步骤：

s41、分别求取步骤s3中产生的每组montecarlo仿真样本的响应值组g1，g2，g3，g4的经验分布函数其中i＝1，2，3，4；其中，n为montecarlo仿真样本的个数，x为经验分布自变量，gij为第j个montecarlo仿真样本带入第i个极限状态方程所求得的响应值，1gij≤x为事件gij≤x的指示函数，当gij≤x时1gij≤x＝1，否则1gij≤x＝0；

s42、将步骤s3中每组样本点响应值组g1，g2，g3，g4带入各自的经验分布函数中，得到新的响应值组y1，y2，y3，y4，其中i＝1，2，3，4；

s43、两两绘制出步骤s42所得新的响应值组y1，y2，…，ym的散点图，观察分析每两种失效模式之间的相关关系，选出备选的copula函数类型：当散点图图像呈现椭圆分布的形状时，采用椭圆族copula函数(正态copula、t-copula等)；当散点图图像呈现上下尾相关的形状时，采用阿基米德族copula函数(gumbelcopula、claytoncopula和frankcopula等)；若散点图中点十分分散均匀的分布在图中的每处，则说明两种失效模式相关关系很小，不需采用copula函数进行分析，此时视为相互独立即可；

在本实施例中，分别绘制出(y1，y2)、(y1，y3)、(y1，y4)、(y2，y3)、(y2，y4)、(y3，y4)的散点图，如图5至图10所示。由图5至图10可知，电容的电参数(容量与等效串联电阻)恶化失效、igbt模块铝键合线与镀通孔疲劳断裂失效，这两对失效模式具有明显的相关关系，且上下尾逐渐变锐呈现出上下尾相关关系，其他对失效模式并不具有明显的相关关系。所以，选取常用的具有尾部相依特征的二元copula函数：gumbelcopula、claytoncopula和frankcopula作为备选copula函数来描述相关性。

s5、根据步骤s3中产生的montecarlo仿真样本，从步骤s4给出的备选copula函数类型中选择最佳copula函数并确定函数的参数值；本步骤中使用kendall秩相关系数确定copula函数的参数值，并选择平方欧氏距离最小的copula函数作为最优函数，包括以下子步骤：

s51、两两求取步骤s42所得的y1，y2，…，ym的kendall秩相关系数：

其中，i与j分别表示第i种和第j种失效模式；n表示montecarlo仿真样本点数；p与q分别表示第p个和第q个样本点；sign为符号函数，即函数内变量为正时，结果为sign＝1，为负时，结果为sign＝-1；gip即表示第i种失效模式的第p个样本点的极限状态方程值。

在本实施例中，可以求得由步骤s43得出的具有相关关系的两对主要失效模式的kendall秩相关系数，结果如表1所示。

表1主要失效模式kendall秩相关系数表

s52、通过步骤s51中所得的kendall秩相关系数，使用解析法分别求出每两种失效模式组的备选copula函数的参数值；本实施例中不同copula函数与kendall秩相关系数表达式如表2所示。

表2不同copula函数与kendall秩相关系数表达式

将步骤s51中求得的kendall秩相关系数代入表2所示的解析表达式，求解方程得到copula函数参数θ的估计值，如表3所示。

表3主要失效模式copula函数参数表

s53、对每组失效模式分别求出每一种备选copula函数的平方欧氏距离，选择平方欧氏距离最小的copula函数为最优copula函数；将步骤s52求得的copula函数参数θ的估计值代入表2所示copula函数表达式中，并求取与经验copula函数的平方欧氏距离，经验copula函数公式为：

其中，i与j分别表示第i种和第j种失效模式；n表示montecarlo仿真样本点数；p表示第p个montecarlo仿真样本点；即为第i种失效模式的第p个样本点的经验分布函数。

平方欧氏距离的公式为：

其中，i与j分别表示第i种和第j种失效模式；n表示montecarlo仿真样本点数；p表示第p个montecarlo仿真样本点；gip即表示第i种失效模式的第p个样本点的极限状态方程值；表示第i种失效模式与第j种失效模式在第p个样本点的极限状态方程值的经验copula函数值，cij(gip,gjp)即表示他们的copula函数值。

得到的结果如表4所示。

表4主要失效模式copula函数拟合优度检验平方欧氏距离表

可以看出，在本实施例中，对于电容的电参数(容量与等效串联电阻)恶化这一失效模式对而言，使用claytoncopula函数来描述两者之间的相关关系比较好；对于igbt模块铝键合线脱落断裂与镀通孔疲劳断裂这一失效模式对而言，使用frankcopula函数来描述两者之间的相关关系比较好。

s6、根据步骤s5所得出的copula函数，计算考虑多失效模式相关的工业机器人电气驱动器可靠度与失效概率。

本实施例中，由步骤s4得到的主要失效模式相关关系：电容的电参数(容量与等效串联电阻)恶化、igbt模块铝键合线脱落断裂与镀通孔疲劳断裂，这两对主要失效模式之间具有明显的相关关系，可以计算考虑多失效模式相关的工业机器人电气驱动器的可靠度r为：

r＝r12×r34(15)

其中，r12表示考虑电容的电参数(容量与等效串联电阻)恶化失效相关关系时，两者的可靠度；r34表示考虑igbt模块铝键合线脱落断裂与镀通孔疲劳断裂失效相关关系时，两者的可靠度。

根据步骤s3所得的各个失效模式的概率与步骤s5所选出的最优copula函数，可以求得r12与r34：

其中，表示使用claytoncopula函数且参数值为21.1688、变量为pf1,pf2的函数值；表示使用frankcopula函数且参数值为19.5944、变量为pf3,pf4的函数值；

所以考虑多失效模式相关的工业机器人电气驱动器的可靠度r为：

r＝r12×r34＝0.9808×0.9783＝0.9595(17)

那么考虑多失效模式相关的工业机器人电气驱动器的失效概率pf为：

pf＝1-r＝1-0.9595＝0.0405(18)

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。