基于改进多层级嵌套Copula方法的机电产品可靠性评估

苏 春 汪千程

(东南大学机械工程学院, 南京 211189)

传统的系统可靠性评估方法,如可靠性框图(RBD)、故障树分析(FTA)等,建立在零部件故障相互独立的假设基础上[1].风力机、数控机床、工程机械等机电产品结构功能复杂且工作环境恶劣,零部件、子系统的故障存在多种相关性和耦合关系,给系统可靠性评估带来挑战.

Murthy等[2]考虑多部件系统可靠性评估问题,率先提出故障相关性的概念,将故障相关分为2种类型:① 瞬时失效相关,它是指系统中某一部件故障会以一定概率导致其他部件瞬时失效;② 渐进退化相关,它是指系统中某一部件故障会在一定程度上影响其他部件的退化速率.在此基础上,Nakagawa等[3]提出第3类故障相关性,即冲击损伤相关,它是指系统中某一部件故障会在一定范围诱导其他部件的随机损伤;当损伤累积到一定程度时部件将会失效.

Copula函数因具有描述变量之间非线性、非对称相关关系的能力,近年来在系统可靠性评估等领域受到关注[4-5].An等[6]以某1.5 MW风电齿轮箱为研究对象,考虑故障相关性,基于Copula方法建立了齿轮可靠性评估模型.Okafor等[7]基于系统与部件的失效概率关系图,利用阿基米德Copula推导出并联系统多态可靠性评估模型.为解决多元变量之间故障相关及结构相依的差异性问题和增强模型适用性,Joe等[8]以二元Copula函数为基础,首次提出成对Copula结构(pair-Copula constructions,PCCs)的概念.Einolander等[9]利用5种不同的Copula函数,完成电动汽车充电事件的仿真与相关性分析,较好地保留了充电事件固有的多样性和依赖性.机电产品通常由若干个子系统组成,每个子系统又包括多个零部件,产品故障存在明显的层次相关性.基于多变量偶发索赔分析(multivariate contingent claim analysis)和成对Copula结构(PCC)理论,Dalla Valle等[10]提出一种新的违约概率估计方法,并利用蒙特卡洛仿真获得近似解.Li等[11]考虑多部件系统中子系统、零部件的多层级随机相关性,利用嵌套Lévy Copula具有的非对称相关结构的建模能力,提出一种视情维修策略,完成多部件系统的维修优化.

受系统内部和外部相关性变量的共同作用,机电产品子系统内部、子系统之间存在故障相关性和层次相关性,Copula函数具有测度变量之间相关关系的能力,能提高可靠性评估精度.但是,已有的研究工作还存在以下不足[6-11]:① 多元变量服从的Copula函数不尽相同,采用同种或同族Copula函数表征变量之间的相关关系存在较大的局限性,无法准确反映变量之间复杂的依赖关系.② 传统嵌套Copula方法在多层级系统的可靠性评估中,联合分布C(·)与实际的联合概率通常存在一定误差.在分层处理过程中,累积误差会降低可靠性评估的准确性,甚至得到错误的分析结果.

本文通过对边缘分布的逐层更新,提出一种改进多层级嵌套Copula的可靠度评估方法.将机电系统分解为子系统级和零部件级,根据故障率模型拟合各零部件的故障概率并作为其边缘分布;通过计算边缘分布之间的秩相关系数,确定零部件之间的故障相关性;再利用改进多层级嵌套Copula完成边缘分布的逐层更新和参数估计,建立零部件级的故障相依结构,通过递归得到各子系统的可靠性模型.重复上述步骤,建立子系统级的故障相依结构,得到系统的可靠性模型.

1 可靠性界限理论

机电产品的结构和功能日趋复杂,科学的可靠性评估是产品安全有效完成既定工作的重要保证.在工作载荷与环境应力的共同作用下,零部件以及系统性能会逐渐退化,直至发生故障.实际上,系统故障的发生具有显著的层次性和相关性.现有的可靠性评估多基于单一层级的可靠性分析模型,未能考虑复杂机电产品的层次性和相关性特征.

根据产品的结构和功能组成,可以将机电产品分解为系统级、子系统级和零部件级等多个层级.如图1所示,机电产品S由多个子系统S1,S2,…,SD组成,其中D为子系统的数量;每个子系统又包含MD个零部件.例如,子系统S2中包括S2,1,S2,2,…,S2,M2个零部件.

图1 机电产品多层级结构组成示意图

根据系统可靠性理论,机电产品可靠度模型可以抽象为各子系统的串联结构.若零部件可靠度已知,根据贝叶斯网络推理即可求得各子系统和系统的可靠度模型,即

(1)

(2)

式中,RS、RS1分别为机电产品(系统)和子系统S1的可靠度;其余变量符号的含义以此类推.

根据可靠性界限理论,系统可靠度的下限Rmin以独立性假设理论为基础,即不考虑子系统之间的故障相关性;其上限Rmax以最薄弱环节理论为基础,即假设子系统之间的故障完全相关[12-13].目前,机电产品可靠性研究多建立在各部件故障与否相互独立的假设基础上,未考虑零部件、子系统的故障相关性.在工程实际中,子系统之间存在一定的依赖关系.因此,系统实际的可靠度Rr应介于Rmin和Rmax之间,即

Rmin(t)=P(TS1>t,TS2>t,…,TSD>t)=

P(TS1>t)P(TS2>t)…P(TSD>t)=

(3)

Rmax(t)=P(TS1>t,TS2>t,…,TSD>t)=

min[RS1(t),RS2(t),…,RSD(t)]

(4)

Rmin(t)

(5)

2 基于多层级嵌套Copula的可靠性建模

Copula函数能够描述变量之间非线性、非对称的相关关系.将机电产品每个零部件的寿命视作低维边缘分布,子系统的寿命视为每个零部件寿命的联合分布,利用Copula函数可实现子系统可靠性评估.同理,将每个子系统的寿命视作低维边缘分布,系统的寿命视为联合分布,利用Copula函数可以实现系统可靠性评估.

2.1 Copula函数及Sklar定理

设F1(x1),F2(x2),…,FN(xN)分别为随机变量x1,x2,…,xN定义域内连续的一元分布函数,N为随机变量的个数.令un=Fn(xn)为随机变量xn(n=1,2,…,N)对应的边缘分布函数,服从[0,1]的均匀分布,则将具有以下性质的多元分布函数C(u1,u2,…,uN)称为N元Copula函数[14]:

1)C(u1,u2,…,uN)的定义域为[0,1]N.

2)C(u1,u2,…,uN)具有N维递增的零基面.

3)C(u1,u2,…,uN)的边缘分布函数Cn(un)满足

Cn(un)=C(1,…,1,un,1,…,1)=un

(6)

Sklar定理指出,Copula函数可以将多个边缘分布连接成一个联合分布,若F(x1,x2,…,xN)为联合分布函数,则存在唯一的Copula函数且满足[14]

F(x1,x2,…,xN)=C[F1(x1),F2(x2),…,FN(xN)]=

C(u1,u2,…,uN)

(7)

2.2 多层级嵌套Copula的定义

令NL表示第L层边缘分布的个数,N元L(L≥1)层嵌套CopulaCL,NL(CL-1,NL-1,1,CL-1,NL-1,2,…,CL-1,NL-1,NL)定义在[0,1]N→[0,1]上,且具有以下性质:

1)CL,NL(CL-1,NL-1,1,CL-1,NL-1,2,…,CL-1,NL-1,NL)具有NL维递增的零基面.

2)当L>1时,CL,NL(CL-1,NL-1,1,CL-1,NL-1,2,…,CL-1,NL-1,NL)的任意边缘分布CL-1,NL-1,a(a=1,2,…,NL)满足下述条件,且至少存在一个边缘分布CL-1,Nl-1,NL(b∈[1,2,…,NL])满足条件①.

①CL-1,NL-1,a为NL-1,a元L-1层嵌套Copula函数.

②CL-1,NL-1,a为Nl,a元l(1≤l

③CL-1,NL-1,a为连续的一元分布函数u1=F1(x1).

3) 当L=1时,C1,a为N1,a元Copula函数C(u1,u2,…,uN1,a).

多层级嵌套Copula利用递归方法,其各层级Copula均可选取拟合效果最好的任意Copula函数.若所有Copula函数均为阿基米德族Copula,则称多层级嵌套Copula为嵌套阿基米德Copula.

2.3 多层级嵌套Copula的建模

2.3.1 确定备选Copula函数

在可靠性评估中,常用的Copula函数包括椭圆族Copula和阿基米德族Copula.其中,椭圆族Copula可分为Gaussian-Copula和t-Copula等;阿基米德族Copula包括Gumbel-Copula、Clayton-Copula和Frank-Copula等.边缘分布为u1,u2,…,uN的N元Gaussian-Copula分布函数为[4]

CGa[(u1,u2,…,uN);θ]=Φθ[Φ-1(u1),

Φ-1(u2),…,Φ-1(uN)]

(8)

式中,θ为相关系数，用来衡量各个变量之间的相关程度;Φθ表示相关系数为θ的N元正态分布;Φ-1(·)为正态分布Φ(·)的逆函数,θ∈[-1,1].

边缘分布为u1,u2,…,uN的N元t-Copula分布函数为[4]

(9)

阿基米德族Copula对变量的边缘分布类型没有特殊限制,能够灵活地描述各类条件下的相关性问题.边缘分布为u1,u2,…,uN的N元阿基米德族Copula分布函数为[4]

C(u1,u2,…,uN)=φ[φ-1(u1)+

φ-1(u2)+…+φ-1(uN)]

(10)

式中,φ为阿基米德Copula的生成元(generator),不同的生成元φ构成不同的阿基米德Copula;φ-1(·)为φ(·)的逆函数,由其生成元完全确定.常用阿基米德族Copula的形式及其参数取值范围如表1所示.

表1 常用阿基米德族Copula的主要形式和参数取值范围

机电产品零部件数量多,不同零部件所服从的Copula函数不尽相同.为准确反映变量之间的相关关系,本文确定以上述5种Copula函数作为备选模型.

2.3.2 Copula函数的相关性度量

在数理统计中,常用Pearson线性相关系数评价2个变量的相关程度.对于考虑边缘分布的Copula模型,变量之间多为非线性相关关系.Dißmann等[15]提出利用Kendall或Spearman秩相关系数来评价变量之间的相关关系.

设连续型随机变量Y和Z的分布函数分别为G(y)和H(z),Copula函数为C(u,v),其中u=G(y),v=H(z),则Kendall秩相关系数τ可以由Copula函数表示为[9]

(11)

设连续型随机变量Y和Z的分布函数分别为G(y)和H(z),Copula函数为C(u,v),其中u=G(y),v=H(z),则Spearman秩相关系数ρ可以由Copula函数表示为[9]

(12)

秩相关系数与Copula函数的参数具有一定的对应关系,通过计算边缘分布之间的秩相关系数即可确定零部件之间的故障相关性.

2.3.3Copula函数的拟合优度检验

不同的Copula可以描述不同的相依结构,根据拟合优度可以选择最能准确刻画变量之间相关关系的Copula函数.作为联合分布函数,F(x1,x2,…,xN)可以分解为一系列成对Copula(pair-Copula)和边缘分布.其中,每组pair-Copula都可以通过赤池信息准则(Akaike information criteria,AIC)来选择多元变量两两之间最优的Copula函数类型.AIC准则定义式为[16]

AIC=-2log(EML)+2W

(13)

式中,EML为极大似然估计;log(EML)为对数似然函数;W为未知参数的个数.AIC的数值越小,则表示所构建的模型拟合度越高.

2.3.4Copula函数的可靠性建模

Copula函数可以连接任意2个随机变量来构造联合分布,准确刻画变量之间的非线性相关关系.以子系统Sd为例,根据Md个零部件建立的故障概率模型FSd和可靠性模型RSd可表示为[13]

FSd(t)=P(TSd,1≤t,TSd,2≤t,…,TSd,Md≤t)=

(-1)MdC(ud,1,ud,2,…,ud,Md)

(14)

(15)

其中,

E2=Ci1,i2(ud,i1,ud,i2)

E3=Ci1,i2,i3(ud,i1,ud,i2,ud,i3)

⋮

Ek=Ci1,i2,…,ik(ud,i1,ud,i2,…,ud,ik)

根据零部件级确定的故障相依结构和相关参数,递归可得到各子系统的可靠性模型.同理,机电产品(系统)S的故障概率模型FS和可靠度模型RS可分别表示为

(16)

(17)

其中,

Q2=Cj1,j2(uj1,uj2)

Q3=Cj1,j2,j3(uj1,uj2,uj3)

⋮

Qk=Cj1,j2,…,jk′(uj1,uj2,…,ujk′)

3 基于改进多层级嵌套Copula的分析

传统嵌套Copula方法将联合分布C(·)作为后续嵌套的新变量,在处理多变量系统时存在误差,将影响可靠性评估的精度.本文通过对边缘分布的逐层更新,提出一种改进多层级嵌套Copula的可靠性评估方法.

设机电系统有D个子系统,每个子系统有MD个零部件,xd,md表示第d个子系统中第md个零部件(d=1,2,…,D;md=1,2,…,MD)所对应的随机变量.基于改进多层级嵌套Copula模型的系统可靠性建模与求解步骤如下:

① 采集零部件的故障数据和试验数据,分析得到各个零部件对应的故障概率,并作为其边缘分布,令ud,md=Fd,md(xd,md),此时各子系统的边缘分布为{u1,1,…,u1,m1,…,u1,M1},…,{ud,1,…,ud,md,…,ud,Md},…,{uD,1,…,uD,mD,…,uD,MD}.

② 以子系统Sd为例,分析其Md个零部件的故障相关性,计算两两变量之间的Kendall秩相关系数τ,并根据τmax确定相关性最强的ud,ip,ud,iq(1≤ip,iq≤MD且ip≠iq)作为最底层变量.

③ 利用ud,ip,ud,iq两个边缘分布构造联合分布,估计最底层备选Copula的参数值Cd,1(ud,ip,ud,iq;θd,1),并基于AIC准则进行拟合优度检验确定最优Copula函数.

⑤ 利用Fd,1和ud,ik两个边缘分布构造联合分布,估计倒数第2层备选Copula的参数值Cd,2(Fd,1,ud,ik;θd,2),并基于AIC准则进行拟合优度检验确定最优Copula函数.

⑥ 将优选出的倒数第2层Copula记作Cd,2,根据Copula可靠性模型求取实际概率Fd,2,实现倒数第2层变量边缘分布的更新.重复步骤④和步骤⑤,直至组内只有1个Copula函数,此时零部件级故障相依结构构建完毕,共有Md-1层Copula.

⑦ 由步骤⑥递归得到各子系统的故障概率Fd(d=1,2,…,D),以此作为边缘分布.重复步骤②～步骤⑥,直至组内只有1个Copula函数,此时子系统级故障相依结构构建完毕,共有D-1层Copula,递归可得机电系统的故障概率F与可靠度R.

4 案例分析

本节以某型号浮式海上风力机为研究对象,完成故障相关条件下多层级系统的可靠度评估,系统结构如图2所示.

根据风力机结构和功能,重点关注变桨系统、齿轮箱和发电机3个子系统,子系统级和零部件级定义如表2所示[17].

子系统和零部件的寿命分布(或故障概率分布)是系统可靠性评估的基础.考虑到此类装备的故障数据特征,本文采用指数模型描述其故障概率,即

P(t)=1-e-λt

(18)

式中,λ为零部件的故障率;P(t)为零部件在t时刻的故障概率.

图2 浮式海上风力机结构功能图

表2 风力机子系统级和零部件级的定义

根据上述模型计算风力机零部件的故障概率并作为其边缘分布.由于S1,3和S3,4的故障率极低,在可靠性分析时不予考虑.

4.1 风力机零部件级的可靠性建模与求解

首先,通过pair-Copula方法对各子系统零部件两两进行相关性分析,采用5种备选Copula函数求解未知参数,计算零部件之间的Kendall秩相关系数并通过AIC准则进行拟合优度检验,其结果如表3所示.

由表3可知,变桨系统中u1,1-u1,2两变量的相关性为0.622 0,齿轮箱中u2,1-u2,2的相关性最高,发电机中u3,1-u3,2的相关性最高.因此,零部件级故障相依结构最底层变量分别为(u1,1,u1,2)、(u2,1,u2,2)和(u3,1,u3,2).根据AIC判别结果拟合最优t-Copula函数C1,1(u1,1,u1,2)、C2,1(u2,1,u2,2)和C3,1(u3,1,u3,2),并由式(14)计算最底层变量的实际联合概率F1,1、F2,1和F3,1,实现齿轮箱和发电机最底层变量边缘分布的更新.

其次,将F2,1和F3,1作为零部件级倒数第2层Copula函数的新变量,与子系统内剩余变量两两进行相关性分析求解未知参数,其分析结果如表4所示.

表3 零部件级最底层变量分析结果

表4 零部件级倒数第2层变量的分析结果

由表4可知,齿轮箱中F2,1-u2,3两变量的相关性为0.616 9,发电机中F3,1-u3,3的相关性最高,因此零部件级故障相依结构倒数第2层变量分别为(F2,1,u2,3)和(F3,1,u3,3),根据AIC判别结果拟合最优t-Copula函数C2,2(F2,1,u2,3)和C3,2(F3,1,u3,3),并计算倒数第2层变量的实际联合概率F2,2和F3,2,实现发电机倒数第2层变量边缘分布的更新.

最后,将F3,2作为零部件级最顶层Copula函数的新变量,与子系统内剩余变量进行相关性分析来求解未知参数,其分析结果如表5所示.

表5 零部件级最顶层变量分析结果

由表5可知,发电机中F3,2-u3,5两变量的相关性为0.425 6,根据AIC判别结果拟合最优t-Copula函数C3,3(F3,2,u3,5),计算最顶层变量的实际联合概率F3,3.综上得到风力机零部件级故障相依结构示意图,如图3所示.

图3 风力机零部件级故障相依结构示意图

由图3可知,根据零部件级的可靠性建模可以得到子系统的故障概率分别为F1,1、F2,2和F3,3,并且改进多层级嵌套Copula方法仅适用于多层级系统(L≥2)的可靠性建模.

针对变桨系统(L=1)的可靠性评估,本文将多层级嵌套Copula方法与嵌套阿基米德Copula方法计算的可靠度进行对比分析.针对齿轮箱和发电机(L≥2)的可靠性评估,将改进多层级嵌套Copula方法与多层级嵌套Copula方法计算的可靠度进行对比分析.绘制风力机各子系统在独立性假设理论、最薄弱环节理论等评估方法下的可靠度曲线,如图4所示.

根据可靠性界限理论,子系统可靠度应介于独立性假设理论和最薄弱环节理论得到的可靠度之间.由图4(a)可知,针对多变量单层级系统S1,多层级嵌套Copula方法可以选择对数据拟合效果最好的t-Copula进行分析,而嵌套阿基米德Copula方法由于备选模型的局限性只能选择拟合效果次之的Gumbel-Copula,计算得到的可靠度难以准确反映变量之间的故障相关性,会影响可靠度评估的精度.

(a) 风力机变桨

(b) 风力机齿轮箱

(c) 风力机发电机

由图4(b)和(c)可知,针对多变量多层级系统S2和S3,改进多层级嵌套Copula方法和传统多层级嵌套Copula方法均可选择对数据拟合效果最好的t-Copula进行分析.然而,前者在每层Copula函数C(·)的基础上先对变量的边缘分布逐层更新,计算变量的实际联合概率,并将其作为后续嵌套的新变量(见图3),可以消除传统方法直接以C(·)作为后续嵌套新变量存在的误差,提高了故障相关条件下多层级系统的可靠度评估精度.

4.2 风力机子系统级的可靠性建模与求解

在零部件级可靠性模型基础上,通过递归得到子系统的故障概率并作为其边缘分布.重复上述步骤,通过pair-Copula方法对各子系统两两进行相关性分析并求解未知参数,计算子系统之间的Kendall秩相关系数,采用AIC准则完成拟合优度检验,结果如表6所示.

由表6可知,风力机中F1,1-F3,3的相关性最高.因此,子系统级故障相依结构的最底层变量为(F1,1,F3,3).根据AIC判别结果拟合最优t-Copula函数C1(F1,1,F3,3),并由式(16)计算最底层变量的实际联合概率F1,实现风力机最底层变量边缘分布的更新.

同理,将F1作为子系统级最顶层Copula函数的新变量,与系统内剩余变量进行相关性分析并求解未知参数,其分析结果如表7所示.

表6 子系统级最底层变量分析结果

表7 子系统级最顶层变量分析结果

由表7可知,风力机中F1-F2,2的相关性为0.477 6.根据AIC判别结果拟合最优t-Copula函数C2(F1,F2,2),计算最顶层变量的联合故障概率,得到风力机子系统级故障相依结构示意图,如图5所示.

图5 风力机子系统级故障相依结构示意图

由图5可知,根据子系统级的可靠性建模得到系统故障概率F,对比改进多层级嵌套Copula方法与多层级嵌套Copula方法计算的可靠度,绘制风力机在独立性假设理论、最薄弱环节理论等方法下的系统可靠度曲线,如图6所示.

由图6可知,基于风力机子系统级的可靠性建模,利用多层级嵌套Copula方法可以得到风力机可靠度曲线.由于三级(零部件级-子系统级-系统级)分层处理的误差累积,多层级嵌套Copula方法

图6 基于不同方法的风力机可靠度曲线对比

得到的可靠度会超出可靠性界限理论的上限,难以准确评估风力机可靠度.本文提出的改进多层级嵌套Copula方法,通过逐层更新边缘分布能消除多层级嵌套Copula方法在分层处理时存在的误差,所得到的可靠度曲线满足可靠性界限理论,符合系统可靠度模型基本规律,能有效分析故障相关条件下多元变量的层次依赖关系.

5 结论

1) 文中提出一种改进多层级嵌套Copula的机电产品可靠度评估方法,将产品分解为系统级、子系统级和零部件级等3个层级.

2) 通过计算边缘分布之间的秩相关系数,确定零部件之间的故障相关性.利用改进多层级嵌套Copula,逐层更新边缘分布并完成参数估计,通过递归建立系统可靠性评估模型.

3) 除串联关系之外,子系统或零部件之间还存在并联、表决等连接关系,后续研究中可以考虑更多的连接关系.Copula函数种类繁多,后续可以进一步改进和优化Copula函数,以更准确地刻画零部件之间的相关性.此外,机电产品可靠度与环境因素的相关性也值得研究.