CN113033002B - 一种针对复杂维纳退化过程模型参数的无偏估计方法 - Google Patents

一种针对复杂维纳退化过程模型参数的无偏估计方法 Download PDF

Info

Publication number
CN113033002B
CN113033002B CN202110326534.3A CN202110326534A CN113033002B CN 113033002 B CN113033002 B CN 113033002B CN 202110326534 A CN202110326534 A CN 202110326534A CN 113033002 B CN113033002 B CN 113033002B
Authority
CN
China
Prior art keywords
degradation
formula
wiener
unbiased
parameter
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110326534.3A
Other languages
English (en)
Other versions
CN113033002A (zh
Inventor
唐圣金
王凤飞
孙晓艳
司小胜
于传强
叶辉
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Rocket Force University of Engineering of PLA
Original Assignee
Rocket Force University of Engineering of PLA
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Rocket Force University of Engineering of PLA filed Critical Rocket Force University of Engineering of PLA
Priority to CN202110326534.3A priority Critical patent/CN113033002B/zh
Publication of CN113033002A publication Critical patent/CN113033002A/zh
Application granted granted Critical
Publication of CN113033002B publication Critical patent/CN113033002B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/13—Differential equations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/04—Ageing analysis or optimisation against ageing

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • Operations Research (AREA)
  • Evolutionary Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Computing Systems (AREA)
  • Probability & Statistics with Applications (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Complex Calculations (AREA)

Abstract

本发明公开了一种针对复杂维纳退化过程模型参数的无偏估计方法，属于可靠性工程技术领域，该方法包括步骤：A、建立设备复杂性能退化模型；B、估计复杂维纳退化过程模型参数；C、计算无偏参数估计。本发明给出了一种针对复杂维纳退化过程模型参数的无偏估计方法，分别针对非线性维纳过程模型、带测量误差的线性维纳过程模型和不定间隔线性维纳过程模型，给出了直接极大化似然估计方法的部分解析参数估计结果，同时给出了相应参数估计的无偏估计和经验无偏估计。为基于复杂维纳过程的无偏估计提供了有力的理论依据，提高了设备剩余寿命预测的精度，从而减少相应的维修经费开支，避免不恰当维修造成的经济损失，具有较好的工程应用价值。

Description

一种针对复杂维纳退化过程模型参数的无偏估计方法
技术领域
本发明属于可靠性工程技术领域，主要内容为一种针对复杂维纳退化过程模型参数的无偏估计方法。
背景技术
维纳过程是现有剩余寿命预测研究中重要的随机模型，常见的复杂维纳过程模型有三种：非线性维纳过程、带测量误差的线性维纳过程和不定间隔线性维纳过程。基于直接极大化似然估计一直没有得到复杂维纳过程模型的解析参数估计结果，以至于复杂维纳过程相关的无偏估计进展缓慢。
发明内容
本发明的目的是提供一种针对复杂维纳退化过程模型参数的无偏估计方法，得到复杂维纳过程模型的解析解和无偏估计，以提高寿命预测精度的方法。
本发明采用如下技术方案来实现的：
一种针对复杂维纳退化过程模型参数的无偏估计方法，包括以下步骤：
1)建立设备复杂性能退化模型
针对不同的退化过程分别建立相应的复杂维纳退化过程模型，得到需要求解的未知参数，即非线性参数θ，扩散系数σB，测量误差
Figure BDA0002994873130000011
和漂移系数的先验信息
Figure BDA0002994873130000012
2)估计复杂维纳退化过程模型参数
利用直接极大化似然估计方法求解各模型中未知参数的解析解，便于进一步分析参数估计的性质；
3)计算无偏参数估计
基于上述未知参数的解析解，分别计算各参数估计对应的期望，在此基础上，提出相应的无偏参数估计，得到最优参数估计结果。
本发明进一步的改进在于，步骤1)的具体实现方法如下：
1)基于非线性维纳过程的退化过程Y(t)表示如下：
Y(t)＝y0+λΛ(t；θ)+σBB(t) (1)
其中，y0表示初始退化状态；λ为漂移系数，表征设备的退化速率；Λ(t；θ)为关于非线性参数θ的单调连续非线性函数，表征设备退化过程的非线性；σB为扩散系数；B(t)为标准布朗运动，表征设备退化过程的动态特性，不失一般性，令y0＝0，为了表示设备个体之间的单元间异质性，漂移系数λ被认为是随机变量，服从期望为μλ和方差为
Figure BDA0002994873130000021
的正态分布，所以，模型中的未知参数为
Figure BDA0002994873130000022
2)基于带测量误差的线性维纳过程的退化过程Z(t)表示如下：
Figure BDA0002994873130000023
其中，z(0)表示初始退化状态，
Figure BDA0002994873130000024
表示测量误差，通常认为ε独立同分布，且与λ相互独立，所以，模型中的未知参数为
Figure BDA0002994873130000025
3)基于不定间隔线性维纳过程的退化过程X(t)表示如下：
X(t)＝x0+λt+σBB(t) (3)
其中，x0表示初始退化状态，模型中的未知参数为
Figure BDA0002994873130000026
本发明进一步的改进在于，步骤2)中，包括非线性维纳过程、带测量误差的线性维纳过程以及不定间隔线性维纳过程。
本发明进一步的改进在于，非线性维纳过程
假设共有n台设备用于测试，每个设备都在相同的监测点t1,t2,...,tm进行监测，则第i台设备在tj时刻监测到的退化数据yi，j表示如下：
yi，j＝Yi(tj)＝λiΛ(tj；θ)+σBB(tj)，1≤i≤n，1≤j≤m (4)
其中，λi为第i设备的退化速率；令Δyi，j＝Yi(tj)-Yi(tj-1)，那么Δyi，j表示为：
Δyi，j＝Yi(tj)-Yi(tj-1)＝λiΔvj+σBB(Δtj)1≤i≤n，1≤j≤m (5)
其中，Δtj＝tj-tj-1，vj＝Λ(tj；θ)，Δvj＝vj-νj-1＝Λ(tj；θ)-Λ(tj-1；θ)；
令Δv＝(Δv1，Δv2，…，Δvm)′，Δyi＝(Δyi，1，Δyi，2，…，Δyi，m)′，Δt＝(Δt1，Δt2，…，Δtm)′和Y＝(Δy′1，Δy′2，…，Δy′n)′，则Δyi服从均值μλΔv和方差
Figure BDA0002994873130000031
的多元正态分布，其中
Figure BDA0002994873130000032
首先，给出式(7)和式(8)
Figure BDA0002994873130000033
Figure BDA0002994873130000034
其中，
Figure BDA0002994873130000035
那么未知参数
Figure BDA0002994873130000036
的对数似然函数ln L(Θ|Y)为
Figure BDA0002994873130000041
对式(9)先求
Figure BDA0002994873130000042
的偏导，然后令偏导为零，得
Figure BDA0002994873130000043
的受限估计
Figure BDA0002994873130000044
Figure BDA0002994873130000045
将式(11)代入到式(9)，得
Figure BDA0002994873130000046
的轮廓似然函数
Figure BDA0002994873130000047
Figure BDA0002994873130000048
最大化式(12)，得
Figure BDA0002994873130000049
的受限估计
Figure BDA00029948731300000410
其中，
Figure BDA00029948731300000411
所以，
Figure BDA0002994873130000051
的受限估计写为
Figure BDA0002994873130000052
将式(13)代入式(12)，得θ的轮廓似然函数ln L(θ|Y)
Figure BDA0002994873130000053
利用MATLAB的FMINSEARCH函数最大化式(15)，即可得到
Figure BDA0002994873130000054
将其代入式(10)、式(13)和式(14)，最终得到
Figure BDA0002994873130000055
带测量误差的线性维纳过程
假设共有n台设备用于测试，每个设备都在相同的监测点t1，t2，...，tm进行监测，则第i台设备在tj时刻监测到的退化数据zi，j表示如下：
zi，j＝Zi(tj)＝λitj+σBB(tj)+εi，j 1≤i≤n，1≤j≤m (16)
其中，λi为第i设备的退化速率；令Δzi，j＝Zi(tj)-Zi(tj-1)，那么Δzi，j表示为：
Δzi，j＝Zi(tj)-Zi(tj-1)＝λiΔtj+σBB(Δtj)+εi，j-εi，j-1 1≤i≤n，1≤j≤m (17)
其中，Δtj＝tj-tj-1；
令Δzi＝(Δzi，1，Δzi，2，...，Δzi，m)′，Δt＝(Δt1，Δt2，…，Δtm)′和Z＝(Δz′1，Δz′2，...，Δz′n)′，则Δzi服从均值μλΔt和方差
Figure BDA0002994873130000056
的多元正态分布，其中，Ω＝diag(Δtj)，F是由fi，j组成的矩阵
Figure BDA0002994873130000057
令D-Ω+φF，其中，
Figure BDA0002994873130000058
那么，
Figure BDA0002994873130000059
则
Figure BDA00029948731300000510
Figure BDA0002994873130000061
那么
Figure BDA0002994873130000062
的对数似然函数ln L(Θ|Z)为
Figure BDA0002994873130000063
最大化式(21)，得
Figure BDA0002994873130000064
Figure BDA0002994873130000065
Figure BDA0002994873130000066
其中，
Figure BDA0002994873130000067
同样的，
Figure BDA0002994873130000068
利用MATLAB的FMINSEARCH函数得到；
不定间隔线性维纳过程
假设共有n台设备用于测试，第i台设备在监测点
Figure BDA0002994873130000069
进行监测，则第i台设备在ti，j时刻监测到的退化数据xi，j表示如下：
xi，j＝Xi(ti，j)＝λiti，j+σBB(ti，j) 1≤i≤n，1≤j≤mi (25)其中，λi为第i设备的退化速率，令Δxi，j＝Xi(ti，j)-Xi(ti，j-1)，那么Δxi，j表示为：
Δxi，j＝Xi(ti，j)-Xi(ti，j-1)＝λiΔti，j+σB(Δti，j)1≤i≤n，1≤j≤mi (26)
其中，Δti，j＝ti，j-ti，j-1；
令
Figure BDA00029948731300000711
和X＝(Δx′1，Δx′2，...，Δx′n)′，则Δxi服从均值μλΔti和方差
Figure BDA0002994873130000071
的多元正态分布，其中
Figure BDA0002994873130000072
则
Figure BDA0002994873130000073
Figure BDA0002994873130000074
那么
Figure BDA0002994873130000075
的对数似然函数为
Figure BDA0002994873130000076
最大化式(30)，得
Figure BDA0002994873130000077
将式(31)代入到式(30)，得
Figure BDA0002994873130000078
的轮廓似然函数：
Figure BDA0002994873130000079
同样的，
Figure BDA00029948731300000710
利用MATLAB的FMINSEARCH函数最大化式(32)得到。
本发明进一步的改进在于，步骤3)中，包括非线性维纳过程、带测量误差的线性维纳过程以及不定间隔线性维纳过程。
本发明进一步的改进在于，非线性维纳过程
由式(1)可得
Figure BDA0002994873130000081
和
Figure BDA0002994873130000082
且(Δv′Ω-1Δyi)′＝Δv′Ω-1Δyi＝Δy′iΩ-1Δv，那么
Figure BDA0002994873130000083
所以
Figure BDA0002994873130000084
Figure BDA0002994873130000085
另外
Figure BDA0002994873130000086
最后，求得
Figure BDA0002994873130000087
的期望
Figure BDA0002994873130000088
基于式(34)、式(35)和式(36)，给出非线性维纳过程参数的无偏估计：
Figure BDA0002994873130000091
Figure BDA0002994873130000092
Figure BDA0002994873130000093
带测量误差的线性维纳过程
同样，因为
Figure BDA0002994873130000094
所以
Figure BDA0002994873130000095
Figure BDA0002994873130000096
Figure BDA0002994873130000097
又因为
Figure BDA0002994873130000098
所以
Figure BDA0002994873130000099
同样地，给出带测量误差的线性维纳过程参数的无偏估计：
Figure BDA0002994873130000101
Figure BDA0002994873130000102
Figure BDA0002994873130000103
不定间隔线性维纳过程
由于无法计算出参数的解析解，所以这里给出经验无偏估计：
Figure BDA0002994873130000104
Figure BDA0002994873130000105
Figure BDA0002994873130000106
其中，
Figure BDA0002994873130000107
本发明至少具有如下有益的技术效果：
本发明给出了针对复杂维纳过程模型参数的无偏估计方法。首先推导出三种复杂维纳过程模型参数估计的解析解，并给出相应的无偏参数估计，完善了复杂维纳过程的参数估计理论体系。另外，无偏参数估计具有更高的参数估计精度，有效提高了设备剩余寿命预测的准确性。
进一步，分别基于非线性维纳过程、带测量误差的线性维纳过程和不定间隔线性维纳过程建立设备性能退化模型，给出前两种模型参数估计的解析解和后一种模型的求解方法，在此解析解的基础上给出了参数估计的性质，并给出相应的无偏参数估计和经验无偏参数估计。复杂维纳过程参数估计的解析解为后续研究提供了有力的技术理论基础，无偏参数估计也为参数估计提供了一个新的方向，有很好的工程应用价值。
附图说明
图1为激光器的退化数据示意图。
图2为非线性维纳过程的仿真退化数据示意图。
图3为带测量误差的线性维纳过程仿真退化数据示意图。
具体实施方式
以下结合附图和实施例对本发明做出进一步的说明。
基于激光器退化数据和仿真数据来验证本发明的有效性。激光器的性能特征通常通过工作电流表征，当工作电流增加的百分比达到某一设定阈值时，激光器失效。Meeker和Escobar给出了15组某激光器的退化数据，测量时间间隔为250h，每组数据包括16个测量点，具体的退化曲线如图1所示，取该激光器的失效阈值w为10。另外，为了验证基于非线性维纳过程和基于带测量误差的线性维纳过程的无偏参数估计的有效性，给出了相应的仿真退化数据，其退化曲线分别如图2，图3所示。
针对复杂维纳退化过程模型参数的无偏估计方法，它包括以下步骤：
A、建立设备复杂性能退化模型
1)基于非线性维纳过程的退化过程Y(t)表示如下：
Y(t)＝y0+λΛ(t；θ)+σBB(t) (1)
其中，y0表示初始退化状态；λ为漂移系数，表征设备的退化速率；Λ(t；θ)为关于非线性参数θ的单调连续非线性函数，表征设备退化过程的非线性；σB为扩散系数；B(t)为标准布朗运动，表征设备退化过程的动态特性。不失一般性，令y0＝0。为了表示设备个体之间的单元间异质性，漂移系数λ被认为是随机变量，服从期望为μλ和方差为
Figure BDA0002994873130000121
的正态分布。所以，模型中的未知参数为
Figure BDA0002994873130000122
2)基于带测量误差的线性维纳过程的退化过程Z(t)表示如下：
Figure BDA0002994873130000123
其中，z(0)表示初始退化状态，
Figure BDA0002994873130000124
表示测量误差，通常认为ε独立同分布，且与λ相互独立。所以，模型中的未知参数为
Figure BDA0002994873130000125
3)基于不定间隔线性维纳过程的退化过程X(t)表示如下：
X(t)＝x0+λt+σBB(t) (3)
其中，x0表示初始退化状态。模型中的未知参数为
Figure BDA0002994873130000126
B、估计复杂维纳退化过程模型参数
1)非线性维纳过程
假设共有n台设备用于测试，每个设备都在相同的监测点t1，t2，...，tm进行监测。则第i台设备在tj时刻监测到的退化数据yi，j表示如下：
yi，j＝Yi(tj)＝λiΛ(tj；θ)+σBB(tj)，1≤i≤n，1≤j≤m (4)
其中，λi为第i设备的退化速率。令Δyi，j＝Yi(tj)-Yi(tj-1)，那么Δyi，j表示为：
Δyi，j＝Yi(tj)-Yi(tj-1)＝λiΔvj+σBB(Δtj) 1≤i≤n，1≤j≤m (5)
其中，Δtj＝tj-tj-1，vj＝Λ(tj；θ)，Δvj＝vj-νj-1＝Λ(tj；θ)-Λ(tj-1；θ)。
令Δv＝(Δv1，Δv2，…，Δvm)′，Δyi＝(Δyi，1，Δyi，2，…，Δyi，m)′，Δt＝(Δt1，Δt2，…，Δtm)′和Y＝(Δy′1，Δy′2，…，Δy′n)′。则Δyi服从均值μλΔv和方差
Figure BDA0002994873130000127
的多元正态分布，其中
Figure BDA0002994873130000128
首先，给出式(7)和式(8)
Figure BDA0002994873130000131
Figure BDA0002994873130000132
其中，
Figure BDA0002994873130000133
那么未知
Figure BDA0002994873130000134
的对数似然函数ln L(Θ|Y)为
Figure BDA0002994873130000135
对式(9)先求
Figure BDA0002994873130000136
的偏导，然后令偏导为零，得
Figure BDA0002994873130000137
的受限估计
Figure BDA0002994873130000138
Figure BDA0002994873130000139
将式(11)代入到式(9)，得
Figure BDA00029948731300001310
的轮廓似然函数
Figure BDA00029948731300001311
Figure BDA0002994873130000141
最大化式(12)，得
Figure BDA0002994873130000142
的受限估计
Figure BDA0002994873130000143
其中，
Figure BDA0002994873130000144
所以，
Figure BDA0002994873130000145
的受限估计写为
Figure BDA0002994873130000146
将式(13)代入式(12)，得θ的轮廓似然函数ln L(θ|Y)
Figure BDA0002994873130000147
利用MATLAB的FMINSEARCH函数最大化式(15)，即可得到
Figure BDA0002994873130000148
将其代入式(10)、式(13)和式(14)，最终得到
Figure BDA0002994873130000149
2)带测量误差的线性维纳过程
假设共有n台设备用于测试，每个设备都在相同的监测点t1，t2，...，tm进行监测。则第i台设备在tj时刻监测到的退化数据zi，j表示如下：
zi，j＝Zi(tj)＝λitj+σBB(tj)+εi，j 1≤i≤n，1≤j≤m (16)
其中，λi为第i设备的退化速率。令Δzi，j＝Zi(tj)-Zi(tj-1)，那么Δzi，j表示为：
Δzi，j＝Zi(tj)-Zi(tj-1)＝λiΔtj+σBB(Δtj)+εi，j-εi，j-1 1≤i≤n，1≤j≤m (17)
其中，Δtj＝tj-tj-1。
令Δzi＝(Δzi，1，Δzi，2，...，Δzi，m)′，Δt＝(Δt1，Δt2，…，Δtm)′和Z＝(Δz′1，Δz′2，...，Δz′n)′。则Δzi服从均值μλΔt和方差
Figure BDA0002994873130000151
的多元正态分布，其中，Ω＝diag(Δtj)，F是由fi，j组成的矩阵
Figure BDA0002994873130000152
令D＝Ω+φF，其中，
Figure BDA0002994873130000153
那么，
Figure BDA0002994873130000154
则
Figure BDA0002994873130000155
Figure BDA0002994873130000156
那么
Figure BDA0002994873130000157
的对数似然函数ln L(Θ|Z)为
Figure BDA0002994873130000158
最大化式(21)，得
Figure BDA0002994873130000159
Figure BDA0002994873130000161
Figure BDA0002994873130000162
其中，
Figure BDA0002994873130000163
同样的，
Figure BDA0002994873130000164
利用MATLAB的FMINSEARCH函数得到。
3)不定间隔线性维纳过程
假设共有n台设备用于测试，第i台设备在监测点
Figure BDA0002994873130000165
进行监测。则第i台设备在ti，j时刻监测到的退化数据xi，j表示如下：
xi，j＝Xi(ti，j)＝λiti，j+σBB(ti，j) 1≤i≤n，1≤j≤mi (25)
其中，λi为第i设备的退化速率。令Δxi，j＝Xi(ti，j)-Xi(ti，j-1)，那么Δxi，j表示为：
Δxi，j＝Xi(ti，j)-Xi(ti，j-1)＝λiΔti，j+σB(Δti，j) 1≤i≤n，1≤j≤mi (26)
其中，Δti，j＝ti，j-ti，j-1。
令
Figure BDA0002994873130000166
和X＝(Δx′1，Δx′2，…，，Δx′n)′。则Δxi服从均值μλΔti和方差
Figure BDA0002994873130000167
的多元正态分布，其中
Figure BDA0002994873130000168
则
Figure BDA0002994873130000169
Figure BDA00029948731300001610
那么
Figure BDA00029948731300001611
的对数似然函数为
Figure BDA0002994873130000171
最大化式(30)，得
Figure BDA0002994873130000172
将式(31)代入到式(30)，得
Figure BDA0002994873130000173
的轮廓似然函数：
Figure BDA0002994873130000174
同样的，
Figure BDA0002994873130000175
利用MATLAB的FMINSEARCH函数最大化式(32)得到。
C、计算无偏参数估计
1)非线性维纳过程
由式(1)可得
Figure BDA0002994873130000176
和
Figure BDA0002994873130000177
且(Δv′Ω-1Δyi)′＝Δv′Ω-1Δyi＝Δy′iΩ-1Δv，那么
Figure BDA0002994873130000178
所以
Figure BDA0002994873130000179
Figure BDA00029948731300001710
另外
Figure BDA0002994873130000181
最后，求得
Figure BDA0002994873130000182
的期望
Figure BDA0002994873130000183
基于式(34)、式(35)和式(36)，给出非线性维纳过程参数的无偏估计：
Figure BDA0002994873130000184
Figure BDA0002994873130000185
Figure BDA0002994873130000186
2)带测量误差的线性维纳过程
同样，因为
Figure BDA0002994873130000187
所以
Figure BDA0002994873130000188
Figure BDA0002994873130000189
Figure BDA0002994873130000191
又因为
Figure BDA0002994873130000192
所以
Figure BDA0002994873130000193
同样地，给出带测量误差的线性维纳过程参数的无偏估计：
Figure BDA0002994873130000194
Figure BDA0002994873130000195
Figure BDA0002994873130000196
3)不定间隔线性维纳过程
由于无法计算出参数的解析解，所以这里给出经验无偏估计：
Figure BDA0002994873130000197
Figure BDA0002994873130000201
Figure BDA0002994873130000202
其中，
Figure BDA0002994873130000203
基于以上步骤，分别针对非线性维纳过程、带测量误差的线性维纳过程和不定间隔的线性维纳过程，求解模型中的未知参数如表1、表2和表3所示：
表1非线性维纳过程参数估计结果
Figure BDA0002994873130000204
表2带测量误差的线性维纳过程参数估计结果
Figure BDA0002994873130000205
表3不定间隔的线性维纳过程参数估计结果
Figure BDA0002994873130000206
从表1、表2、表3看出，无偏估计值与相比于直接极大化似然估计方法获得估计值更接近于真实值。
综上所述，无偏估计方法优于极大化似然估计方法，本发明提出的方法可以提高模型参数估计的准确性，并提高剩余寿命估计的准确性，进而验证了本发明的有效性。

Claims (4)

1.一种针对复杂维纳退化过程模型参数的无偏估计方法，其特征在于，包括以下步骤：
1)建立设备复杂性能退化模型
针对不同的退化过程分别建立相应的复杂维纳退化过程模型，得到需要求解的未知参数，即非线性参数θ，扩散系数σB，测量误差
Figure FDA0003970768130000011
和漂移系数的先验信息
Figure FDA0003970768130000012
2)估计复杂维纳退化过程模型参数
利用直接极大化似然估计方法求解各模型中未知参数的解析解，便于进一步分析参数估计的性质；
非线性维纳过程
假设共有n台设备用于测试，每个设备都在相同的监测点t1,t2,...,tm进行监测，则第i台设备在tj时刻监测到的退化数据yi,j表示如下：
yi,j＝Yi(tj)＝λiΛ(tj；θ)+σBB(tj),1≤i≤n,1≤j≤m (4)
其中，λi为第i设备的退化速率；令Δyi,j＝Yi(tj)-Yi(tj-1)，那么Δyi,j表示为：
Δyi,j＝Yi(tj)-Yi(tj-1)＝λiΔνj+σBB(Δtj) 1≤i≤n,1≤j≤m (5)
其中，Δtj＝tj-tj-1，νj＝Λ(tj；θ)，Δνj＝νj-νj-1＝Λ(tj；θ)-Λ(tj-1；θ)；
令Δv＝(Δv1,Δv2,…,Δvm)′，Δyi＝(Δyi,1,Δyi,2,…,Δyi,m)′，Δt＝(Δt1,Δt2,…,Δtm)′和Y＝(Δy′1,Δy′2,…,Δy′n)′，则Δyi服从均值μλΔν和方差
Figure FDA0003970768130000013
的多元正态分布，其中
Figure FDA0003970768130000014
首先，给出式(7)和式(8)
Figure FDA0003970768130000015
Figure FDA0003970768130000021
其中，
Figure FDA0003970768130000022
那么未知参数
Figure FDA0003970768130000023
的对数似然函数lnL(Θ|Y)为
Figure FDA0003970768130000024
对式(9)先求μλ,
Figure FDA0003970768130000025
的偏导，然后令偏导为零，得μλ,
Figure FDA0003970768130000026
的受限估计
Figure FDA0003970768130000027
Figure FDA0003970768130000028
将式(11)代入到式(9)，得
Figure FDA0003970768130000029
的轮廓似然函数
Figure FDA00039707681300000210
Figure FDA00039707681300000211
最大化式(12)，得
Figure FDA00039707681300000212
的受限估计
Figure FDA0003970768130000031
其中，
Figure FDA0003970768130000032
所以，
Figure FDA0003970768130000033
的受限估计写为
Figure FDA0003970768130000034
将式(13)代入式(12)，得θ的轮廓似然函数lnL(θ|Y)
Figure FDA0003970768130000035
利用MATLAB的FMINSEARCH函数最大化式(15)，即可得到
Figure FDA0003970768130000036
将其代入式(10)、式(13)和式(14)，最终得到μλ,
Figure FDA0003970768130000037
带测量误差的线性维纳过程
假设共有n台设备用于测试，每个设备都在相同的监测点t1,t2,...,tm进行监测，则第i台设备在tj时刻监测到的退化数据zi,j表示如下：
zi,j＝Zi(tj)＝λitj+σBB(tj)+εi,j 1≤i≤n,1≤j≤m (16)
其中，λi为第i设备的退化速率；令Δzi,j＝Zi(tj)-Zi(tj-1)，那么Δzi,j表示为：
Δzi,j＝Zi(tj)-Zi(tj-1)＝λiΔtj+σBB(Δtj)+εi,j-εi,j-1 1≤i≤n,1≤j≤m (17)
其中，Δtj＝tj-tj-1；
令Δzi＝(Δzi,1,Δzi,2,...,Δzi,m)′，Δt＝(Δt1,Δt2,…,Δtm)′和Z＝(Δz′1,Δz′2,...,Δz′n)′，则Δzi服从均值μλΔt和方差
Figure FDA0003970768130000038
的多元正态分布，其中，Ω＝diag(Δtj)，F是由fi,j组成的矩阵
Figure FDA0003970768130000041
令DΩ＝+φF，其中，
Figure FDA0003970768130000042
那么，
Figure FDA0003970768130000043
则
Figure FDA0003970768130000044
Figure FDA0003970768130000045
那么
Figure FDA0003970768130000046
的对数似然函数lnL(Θ|Z)为
Figure FDA0003970768130000047
最大化式(21)，得
Figure FDA0003970768130000048
Figure FDA0003970768130000049
Figure FDA00039707681300000410
其中，
Figure FDA00039707681300000411
同样的，
Figure FDA00039707681300000412
利用MATLAB的FMINSEARCH函数得到；
不定间隔线性维纳过程
假设共有n台设备用于测试，第i台设备在监测点ti,1,ti,2,…,ti,mi进行监测，则第i台设备在ti,j时刻监测到的退化数据xi,j表示如下：
xi,j＝Xi(ti,j)＝λiti,j+σBB(ti,j) 1≤i≤n,1≤j≤mi (25)
其中，λi为第i设备的退化速率，令Δxi,j＝Xi(ti,j)-Xi(ti,j-1)，那么Δxi,j表示为：
Δxi,j＝Xi(ti,j)-Xi(ti,j-1)＝λiΔti,j+σB(Δti,j) 1≤i≤n,1≤j≤mi (26)
其中，Δti,j＝ti,j-ti,j-1；
令
Figure FDA0003970768130000051
和X＝(Δx′1,Δx′2,…,Δx′n)′，则Δxi服从均值μλΔti和方差
Figure FDA0003970768130000052
的多元正态分布，其中
Figure FDA0003970768130000053
则
Figure FDA0003970768130000054
Figure FDA0003970768130000055
那么
Figure FDA0003970768130000056
的对数似然函数为
Figure FDA0003970768130000057
最大化式(30)，得
Figure FDA0003970768130000058
将式(31)代入到式(30)，得
Figure FDA0003970768130000061
的轮廓似然函数：
Figure FDA0003970768130000062
同样的，
Figure FDA0003970768130000063
利用MATLAB的FMINSEARCH函数最大化式(32)得到；
3)计算无偏参数估计
基于上述未知参数的解析解，分别计算各参数估计对应的期望，在此基础上，提出相应的无偏参数估计，得到最优参数估计结果。
2.根据权利要求1所述的一种针对复杂维纳退化过程模型参数的无偏估计方法，其特征在于，步骤1)的具体实现方法如下：
1)基于非线性维纳过程的退化过程Y(t)表示如下：
Y(t)＝y0+λΛ(t；θ)+σBB(t) (1)
其中，y0表示初始退化状态；λ为漂移系数，表征设备的退化速率；Λ(t；θ)为关于非线性参数θ的单调连续非线性函数，表征设备退化过程的非线性；σB为扩散系数；B(t)为标准布朗运动，表征设备退化过程的动态特性，不失一般性，令y0＝0，为了表示设备个体之间的单元间异质性，漂移系数λ被认为是随机变量，服从期望为μλ和方差为
Figure FDA0003970768130000064
的正态分布，所以，模型中的未知参数为
Figure FDA0003970768130000065
2)基于带测量误差的线性维纳过程的退化过程Z(t)表示如下：
Figure FDA0003970768130000066
其中，z0表示初始退化状态，
Figure FDA0003970768130000067
表示测量误差，ε独立同分布，且与λ相互独立，所以，模型中的未知参数为
Figure FDA0003970768130000068
3)基于不定间隔线性维纳过程的退化过程X(t)表示如下：
X(t)＝x0+λt+σBB(t) (3)
其中，x0表示初始退化状态，模型中的未知参数为
Figure FDA0003970768130000071
3.根据权利要求1所述的一种针对复杂维纳退化过程模型参数的无偏估计方法，其特征在于，步骤3)中，包括非线性维纳过程、带测量误差的线性维纳过程以及不定间隔线性维纳过程。
4.根据权利要求3所述的一种针对复杂维纳退化过程模型参数的无偏估计方法，其特征在于，非线性维纳过程
由式(1)可得
Figure FDA0003970768130000072
和
Figure FDA0003970768130000073
且(Δv′Ω-1Δyi)′＝Δv′Ω-1Δyi＝ΔyiΩ-1Δv，那么
Figure FDA0003970768130000074
所以
Figure FDA0003970768130000075
Figure FDA0003970768130000076
另外
Figure FDA0003970768130000077
最后，求得
Figure FDA0003970768130000078
的期望
Figure FDA0003970768130000081
基于式(34)、式(35)和式(36)，给出非线性维纳过程参数的无偏估计：
Figure FDA0003970768130000082
Figure FDA0003970768130000083
Figure FDA0003970768130000084
带测量误差的线性维纳过程
同样，因为
Figure FDA0003970768130000085
所以
Figure FDA0003970768130000086
Figure FDA0003970768130000087
Figure FDA0003970768130000088
又因为
Figure FDA0003970768130000091
所以
Figure FDA0003970768130000092
同样地，给出带测量误差的线性维纳过程参数的无偏估计：
Figure FDA0003970768130000093
Figure FDA0003970768130000094
Figure FDA0003970768130000095
不定间隔线性维纳过程
由于无法计算出参数的解析解，所以这里给出经验无偏估计：
Figure FDA0003970768130000096
Figure FDA0003970768130000097
Figure FDA0003970768130000098
其中，
Figure FDA0003970768130000099
CN202110326534.3A 2021-03-26 2021-03-26 一种针对复杂维纳退化过程模型参数的无偏估计方法 Active CN113033002B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110326534.3A CN113033002B (zh) 2021-03-26 2021-03-26 一种针对复杂维纳退化过程模型参数的无偏估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110326534.3A CN113033002B (zh) 2021-03-26 2021-03-26 一种针对复杂维纳退化过程模型参数的无偏估计方法

Publications (2)

Publication Number Publication Date
CN113033002A CN113033002A (zh) 2021-06-25
CN113033002B true CN113033002B (zh) 2023-03-21

Family

ID=76474196

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110326534.3A Active CN113033002B (zh) 2021-03-26 2021-03-26 一种针对复杂维纳退化过程模型参数的无偏估计方法

Country Status (1)

Country Link
CN (1) CN113033002B (zh)

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2008032437A (ja) * 2006-07-26 2008-02-14 Toshiba Corp ウェイト算出方法、ウェイト算出装置、アダプティブアレーアンテナ、及びレーダ装置
CN107145645A (zh) * 2017-04-19 2017-09-08 浙江大学 带不确定冲击的非平稳退化过程剩余寿命预测方法
CN107862134A (zh) * 2017-11-06 2018-03-30 河南科技大学 一种考虑自相关测量误差的Wiener过程可靠性分析方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108959676B (zh) * 2017-12-22 2019-09-20 北京航空航天大学 一种考虑有效冲击的退化建模与寿命预测方法
CN109977491B (zh) * 2019-03-06 2020-11-13 北京航空航天大学 一种冲击损伤可恢复条件下的退化建模与寿命预测方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2008032437A (ja) * 2006-07-26 2008-02-14 Toshiba Corp ウェイト算出方法、ウェイト算出装置、アダプティブアレーアンテナ、及びレーダ装置
CN107145645A (zh) * 2017-04-19 2017-09-08 浙江大学 带不确定冲击的非平稳退化过程剩余寿命预测方法
CN107862134A (zh) * 2017-11-06 2018-03-30 河南科技大学 一种考虑自相关测量误差的Wiener过程可靠性分析方法

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
Comparative Study on Simulation of Life and Degradation Tests Based on Wiener Degradation Failure Products;Xin-Yu Ma et al;《IEEE》;20200305;全文 *
The wiener degradation model in reliability analysis;Evgeniya S et al;《IEEE》;20170323;全文 *
基于Wiener过程的DRO加速退化试验方法研究;麻连净等;《装备环境工程》;20100615(第03期);全文 *
基于随机维纳过程的产品加速因子分布确定方法;魏高乐等;《科学技术与工程》;20151018(第29期);全文 *
机载电子设备在线可靠性评估与剩余寿命预测方法研究;王书锋;《中国优秀硕士学位论文全文数据库（电子期刊）工程科技Ⅱ辑》;20150115;第2015年卷(第1期);全文 *
考虑有损伤测量条件下的设备剩余寿命的建模分析;王小林等;《国防科技大学学报》;20111028(第05期);全文 *

Also Published As

Publication number Publication date
CN113033002A (zh) 2021-06-25

Similar Documents

Publication Publication Date Title
CN109829137B (zh) 一种周期应力下非线性退化设备的寿命预测方法及系统
CN111046564B (zh) 两阶段退化产品的剩余寿命预测方法
CN109828220B (zh) 一种锂离子电池健康状态线性评估方法
CN106295003B (zh) 一种基于退化轨迹坐标重构和多元线性回归的锂电池寿命预测方法
CN105372531A (zh) 基于Weibull分布模型的变压器绝缘热老化参数相关性计算方法
CN113033002B (zh) 一种针对复杂维纳退化过程模型参数的无偏估计方法
CN116086537A (zh) 一种设备状态监测方法、装置、设备及存储介质
CN113791351B (zh) 基于迁移学习和差值概率分布的锂电池寿命预测方法
CN108537423B (zh) 一种核电厂对数正态分布可靠性数据处理的Bayes方法
CN111914386A (zh) 一种基于退化模型不确定分析的可靠性评估方法及系统
CN114792053B (zh) 一种基于初值-速率相关退化模型的可靠性评估方法
CN113011036B (zh) 一种针对线性维纳退化过程模型参数的无偏估计方法
CN112748663B (zh) 一种基于数据驱动输出反馈的风电转矩容错控制方法
CN112098782A (zh) 一种基于神经网络的moa绝缘状态检测方法及系统
CN112491038B (zh) 一种配电网设备寿命的概率分布估计方法及系统
CN111044176A (zh) 一种发电机温度异常的监测方法
CN113553232B (zh) 一种通过在线矩阵画像对运维数据进行无监督异常检测的技术
CN109814390B (zh) Miso异因子全格式无模型控制方法
CN114646891B (zh) 一种结合lstm网络和维纳过程的剩余寿命预测方法
CN109782587B (zh) Miso异因子紧格式无模型控制方法
CN104991995A (zh) 纯锡镀层元器件锡须生长失效预测方法与系统
CN117787156A (zh) 一种向量回归与维纳结合的集成电路寿命预测方法
CN116227191A (zh) 结合多项式拟合与维纳过程的锂电池剩余寿命预测方法
CN116029165A (zh) 一种计及雷电影响的电力电缆可靠性分析方法及系统
CN117192415A (zh) 基于误差修正的锂电池剩余寿命自适应预测方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant