CN111597647B - 一种面向工业生产过程的弹簧阻尼系统滤波故障诊断方法 - Google Patents

一种面向工业生产过程的弹簧阻尼系统滤波故障诊断方法 Download PDF

Info

Publication number
CN111597647B
CN111597647B CN202010272258.2A CN202010272258A CN111597647B CN 111597647 B CN111597647 B CN 111597647B CN 202010272258 A CN202010272258 A CN 202010272258A CN 111597647 B CN111597647 B CN 111597647B
Authority
CN
China
Prior art keywords
spring damping
damping system
time
intersection
fault
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
CN202010272258.2A
Other languages
English (en)
Other versions
CN111597647A (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.)
Jiangnan University
Original Assignee
Jiangnan University
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 Jiangnan University filed Critical Jiangnan University
Priority to CN202010272258.2A priority Critical patent/CN111597647B/zh
Publication of CN111597647A publication Critical patent/CN111597647A/zh
Application granted granted Critical
Publication of CN111597647B publication Critical patent/CN111597647B/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/10—Geometric CAD
    • G06F30/17—Mechanical parametric or variational design
    • 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/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Physics & Mathematics (AREA)
  • Geometry (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Testing Of Devices, Machine Parts, Or Other Structures Thereof (AREA)

Abstract

本发明公开了一种面向工业生产过程的弹簧阻尼系统滤波故障诊断方法，属于弹簧阻尼系统故障诊断技术领域。所述方法通过弹簧阻尼系统故障指示信号的数值，确定系统的故障状态，在检测出系统发生故障时，根据测试集合及其
Figure ZY_1
椭球确定弹簧阻尼系统的参数向量θ的具体故障分量，继而根据具体故障分量得到扩展方向，按扩展方向重置支持正多胞体的交集，相对于现有故障检测方法中的全域扩展，本申请所提出的面向工业生产过程的弹簧阻尼系统滤波故障诊断方法收敛速度更快、可以更快的识别出故障参数、从而实时性能更好。

Description

一种面向工业生产过程的弹簧阻尼系统滤波故障诊断方法
技术领域
本发明涉及一种面向工业生产过程的弹簧阻尼系统滤波故障诊断方法，属于弹簧阻尼系统故障诊断技术领域。
背景技术
弹簧阻尼系统是工业生产过程涉及的机械振动系统当中普遍使用的一类结构，在工业生产过程中主要起到吸收和耗散生产过程能量的作用，其吸收和耗散能量的大小，关系到生产过程的安全稳定，所以其系统的安全性将会对整个系统的稳定运行产生重大的影响。一旦弹簧阻尼系统发生故障而未被及时诊断，那么将会影响整个系统的运行。因此，对弹簧阻尼系统进行故障诊断的研究对工业生产过程有着重要的现实意义。
常见的故障诊断方法主要分为基于模型的故障诊断方法和基于数据的故障诊断方法两大类。其中基于数据的故障诊断方法需以大量的数据分析为前提，很难做到实时故障诊断。基于模型的故障诊断方法以系统解析模型为基础，能够实现实时在线的故障诊断；但是一般的基于模型的故障诊断方法中要求系统噪声符合一定概率分布，然而实际应用中，弹簧阻尼系统常受各种外部因素的影响，噪声干扰因素不确定，导致其无法满足这一要求。
集员滤波方法只要求系统噪声有界，对噪声的概率分布没有限定，因此基于集员估计的故障诊断方法能够有效处理噪声不确定系统的故障诊断问题。集员估计方法中是基于可行集的状态来确定系统的故障状态。若检测到可行集为空，则认为系统存在故障；若系统可行集不为空，则认为系统没有故障。集员估计方法中根据近似可行集的空间形状和包围方法的不同，可将其划分为椭球体、全对称多胞体、正多胞体等等。
现有的基于集员滤波理论的故障诊断方法中，在检测出系统发生故障后，若要通过重置参数集合的方法来进一步识别故障，一般是在对参数集合进行全域扩展下进行的，然而全域扩展下的参数集合在参数向量的各个维度都需要再次进行收缩，所以存在着收敛速度慢、故障识别时间长，实时性较差的问题。
发明内容
为了解决目前存在的问题，本发明提供了一种面向工业生产过程的弹簧阻尼系统滤波故障诊断方法，所述方法包括：
根据弹簧阻尼系统的系统方程和信息向量确定弹簧阻尼系统的参数向量对应的椭球集合，进而确定椭球集合对应的支持正多胞体集合和支持正多胞体的交集；
在系统发生故障时，根据测试集合及其
Figure GDA0004135471660000023
椭球确定弹簧阻尼系统的参数向量的具体故障分量继而根据弹簧阻尼系统的参数向量的具体故障分量确定扩展方向，按扩展方向重置支持正多胞体的交集；
根据按扩展方向重置的支持正多胞体的交集，识别弹簧阻尼系统的参数向量中的故障参数。
可选的，所述方法包括：
步骤101，建立弹簧阻尼系统的离散系统模型，确定弹簧阻尼系统的系统方程；
y(k)＝θTφ(k)+e(k)
其中，y(k)为k时刻系统的输出，θ＝[θ1,θ2,θ3,θ4]T为参数向量，k表示离散时间，φ(k)为信息向量，e(k)为弹簧阻尼系统的不确定噪声，且e(k)有界，即|e(k)|≤γ，γ为大于零的常数；
步骤102，获取弹簧阻尼系统在工作状态下的外加控制力和对应的物块的位移，以确定系统的信息向量φ(k)；
步骤103，根据步骤101确定的弹簧阻尼系统的系统方程和步骤102确定的信息向量φ(k)，确定参数向量θ对应的椭球集合和系统故障指示信号f(k)的数值；
步骤104，根据步骤103确定的弹簧阻尼系统参数向量θ对应的椭球集合，确定椭球集合对应的支持正多胞体集合O(k)和支持正多胞体的交集；
步骤105，根据步骤103确定的弹簧阻尼系统的故障指示信号f(k)，确定系统的故障状态和故障时间；
步骤106，若系统发生故障，确定弹簧阻尼系统的参数向量θ的具体故障分量；
步骤107，根据步骤106确定的弹簧阻尼系统的参数向量的具体故障分量确定扩展方向，按扩展方向重置支持正多胞体的交集；
步骤108，根据步骤107中按扩展方向重置的支持正多胞体的交集，识别弹簧阻尼系统的故障参数。
可选的，所述步骤104，根据步骤103确定的弹簧阻尼系统参数向量θ对应的椭球集合，确定椭球集合对应的支持正多胞体集合O(k)和支持正多胞体的交集，包括：
Figure GDA0004135471660000021
Figure GDA0004135471660000022
Figure GDA0004135471660000031
其中，
Figure GDA0004135471660000032
表示k时刻支持正多胞体集合O(k)的参数的上界，
Figure GDA0004135471660000033
表示k时刻支持正多胞体集合O(k)的参数的下界，
Figure GDA0004135471660000034
u∈{1,2,…,n}，<·>为内积函数，n为参数向量θ的维数，k时刻表示第k个离散时刻；
根据支持正多胞体集合O(k)，确定k时刻支持正多胞体的交集：
X(k)＝O(1)∩…∩O(k)＝X(k-1)∩O(k)
即
Figure GDA0004135471660000035
Figure GDA0004135471660000036
Figure GDA0004135471660000037
其中，
Figure GDA0004135471660000038
表示k时刻支持正多胞体的交集X(k)的参数的上界，
Figure GDA0004135471660000039
表示k时刻支持正多胞体的交集X(k)的参数的下界。
可选的，所述步骤106，若系统发生故障，确定弹簧阻尼系统的参数向量θ的具体故障分量，包括：
对k-1时刻支持正多胞体的交集X(k-1)进行n-1维扩展，获得测试集合
Figure GDA00041354716600000310
i∈{1,…,n}，即对于第i个测试集合
Figure GDA00041354716600000311
除掉第i维，都进行扩展；
根据所述k-1时刻的测试集合
Figure GDA00041354716600000312
计算k-1时刻测试集合对应的
Figure GDA00041354716600000313
椭球
Figure GDA00041354716600000314
利用弹簧阻尼系统的系统方程和信息向量，更新k-1时刻的
Figure GDA00041354716600000325
椭球
Figure GDA00041354716600000315
得到k时刻的
Figure GDA00041354716600000316
椭球
Figure GDA00041354716600000317
计算所述k+L时刻
Figure GDA00041354716600000318
椭球
Figure GDA00041354716600000319
空集指示信号
Figure GDA00041354716600000320
根据所述k+L时刻
Figure GDA00041354716600000321
椭球
Figure GDA00041354716600000322
空集指示信号
Figure GDA00041354716600000323
确定弹簧阻尼系统的参数向量θ的具体故障分量。
可选的，所述步骤107，根据步骤106确定的弹簧阻尼系统的参数向量的具体故障分量确定扩展方向，按扩展方向重置支持正多胞体的交集，包括：
若弹簧阻尼系统的参数向量θ中的θi发生故障，按下式更新k-1时刻重置的支持正多胞体的交集Xr(k-1)的第i维参数的上界：
Figure GDA00041354716600000324
若弹簧阻尼系统的参数向量θ中的θi发生故障，按下式更新k-1时刻重置的支持正多胞体的交集Xr(k-1)的第i维参数的下界：
Figure GDA0004135471660000041
若弹簧阻尼系统的参数向量θ中的θi未发生故障，按下式更新k-1时刻重置的支持正多胞体的交集Xr(k-1)的第i维参数的上界：
Figure GDA0004135471660000042
若弹簧阻尼系统的参数向量θ中的θi未发生故障，按下式更新k-1时刻重置的支持正多胞体的交集Xr(k-1)的第i维参数的下界：
Figure GDA0004135471660000043
根据所述k-1时刻支持正多胞体的交集Xr(k-1)每维参数的上界和下界，按下式获得k-1时刻按扩展方向重置的支持正多胞体的交集Xr(k-1)：
Figure GDA0004135471660000044
其中，
Figure GDA0004135471660000045
为第i维参数θi的最大变化范围，i∈{1,…,n}，r为上标，带有此上标的参数和集合表示重置后的参数和集合。
可选的，所述根据步骤107中按扩展方向重置的支持正多胞体的交集，识别弹簧阻尼系统的故障参数，包括：
根据所述k-1时刻按扩展方向重置的支持正多胞体的交集Xr(k-1)，计算k-1时刻按扩展方向重置的支持正多胞体的交集Xr(k-1)对应的
Figure GDA0004135471660000046
椭球Er(k-1)；
利用k时刻物块的位移y(k)、信息向量φ(k)，k-1时刻的
Figure GDA0004135471660000047
椭球Er(k-1)的中心θcr(k-1)、轴信息矩阵Pr(k-1)和系统不确定噪声的边界γ，重新更新k时刻弹簧阻尼系统参数向量对应的椭球集合E(k)的中心θc(k)和轴信息矩阵P(k)；
根据所述重新更新的椭球集合E(k)，重新更新k时刻椭球集合E(k)对应的支持正多胞体集合O(k)和支持正多胞体的交集X(k)；
根据所述支持正多胞体的交集X(k)每一维的上界和下界，确定弹簧阻尼系统的故障参数向量
Figure GDA0004135471660000048
Figure GDA0004135471660000049
其中
Figure GDA0004135471660000051
Figure GDA0004135471660000052
为k时刻支持正多胞体的交集X(k)的第i维对应的参数的上界，
Figure GDA0004135471660000053
为k时刻支持正多胞体的交集X(k)的第i维对应的参数的下界。
可选的，所述步骤103，根据步骤101确定的弹簧阻尼系统的系统方程和步骤102确定的信息向量φ(k)，确定参数向量θ对应的椭球集合和系统故障指示信号f(k)的数值，包括：
按下面两式确定k时刻仿射变换后
Figure GDA0004135471660000054
正交的两个平行超平面的
Figure GDA0004135471660000055
坐标，即α(k)和
Figure GDA0004135471660000056
其中
Figure GDA0004135471660000057
为仿射变换后的φ：
Figure GDA0004135471660000058
Figure GDA0004135471660000059
若α(k)≥1或
Figure GDA00041354716600000510
则故障指示信号f(k)＝1，表明弹簧阻尼系统在k时刻发生故障；
θc(k)＝θc(k-1)，
P(k)＝P(k-1)；
若α(k)≤1且
Figure GDA00041354716600000511
则故障指示信号f(k)＝0，表明弹簧阻尼系统在k时刻未发生故障。
可选的，若
Figure GDA00041354716600000512
按下面两式更新k时刻椭球集合E(k)的中心θc(k)和轴信息矩阵P(k)：
θc(k)＝θc(k-1)，
P(k)＝P(k-1)；
若两者还同时满足
Figure GDA00041354716600000513
则
Figure GDA00041354716600000514
在上述条件下，若|μ(k)|＞ρ，则
Figure GDA00041354716600000515
Figure GDA00041354716600000516
Figure GDA0004135471660000061
Figure GDA0004135471660000062
Figure GDA0004135471660000063
在上述条件下，若|μ(k)|≤ρ，则
Figure GDA0004135471660000064
τ(k)＝0，
σ(k)＝nα2，
Figure GDA0004135471660000065
之后按下面两式更新k时刻椭球集合E(k)的中心θc(k)和轴信息矩阵P(k)
Figure GDA0004135471660000066
Figure GDA0004135471660000067
其中n为参数向量θ的维数，μ(k)为α(k)和
Figure GDA0004135471660000068
的平均值，ρ为大于零的数，设定ρ＝10-6，τ(k)为k时刻仿射变换后
Figure GDA0004135471660000069
沿
Figure GDA00041354716600000610
的中心坐标，σ(k)为
Figure GDA00041354716600000611
沿
Figure GDA00041354716600000612
半轴的平方，
Figure GDA00041354716600000613
为
Figure GDA00041354716600000614
正交于
Figure GDA00041354716600000615
的半轴的平方，ε(k)、b(k)、α为中间变量，
Figure GDA00041354716600000616
为仿射变换后的E，
Figure GDA00041354716600000617
为仿射变换后的φ。
可选的，所述步骤102，获取弹簧阻尼系统在工作状态下的外加控制力和对应的物块的位移，以确定系统的信息向量φ(k)，包括：
在预定时间范围内，获取弹簧阻尼系统在工作状态下的外加控制力和对应的物块的位移；
将所得的外加控制力和对应的物块的位移的数据代入下式中：
φ(k)＝[-y(k-1),-y(k-2),u(k-1),u(k-2)]T
确定出弹簧阻尼系统的信息向量φ(k)；k的取值范围为1至N，k为整数。
本申请的另一个方面还提供一种弹簧阻尼系统滤波故障诊断系统，所述系统采用上述面向工业生产过程的弹簧阻尼系统滤波故障诊断方法对弹簧阻尼系统进行故障诊断。
本发明有益效果是：
本申请通过弹簧阻尼系统故障指示信号的数值，确定系统的故障状态，在检测出系统发生故障时，根据测试集合及其
Figure GDA0004135471660000071
椭球确定弹簧阻尼系统的参数向量θ的具体故障分量，继而根据具体故障分量得到扩展方向，按扩展方向重置支持正多胞体的交集，相对于现有故障检测方法中的全域扩展，本申请所提出的面向工业生产过程的弹簧阻尼系统滤波故障诊断方法收敛速度更快、可以更快的识别出故障参数、从而实时性能更好。
附图说明
为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1是本发明一个实施例中公开的一种面向工业生产过程的弹簧阻尼系统滤波故障诊断方法流程图。
图2是一种弹簧系统结构图。
图3是本发明一个实施例中公开的两种算法下的一种弹簧系统的参数分量θ1的上下界变化曲线图。
图4是本发明一个实施例中公开的两种算法下的一种弹簧系统的参数分量θ2的上下界变化曲线图。
图5是本发明一个实施例中公开的两种算法下的一种弹簧系统的参数分量θ3的上下界变化曲线图。
图6是本发明一个实施例中公开的两种算法下的一种弹簧系统的参数分量θ4的上下界变化曲线图。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地详细描述。
实施例一：
本实施例提供一种面向工业生产过程的弹簧阻尼系统滤波故障诊断方法，参见图1，所述方法包括：
步骤101，建立弹簧阻尼系统的离散系统模型，确定弹簧阻尼系统的系统方程。
弹簧阻尼系统的结构图如图2所示，其中M为物块的质量，C为摩擦系数，K是弹簧的弹性系数，u(t)为t时刻的外加控制力，是系统的输入，y(t)为t时刻的物块的位移，是系统的输出。
根据此模型结构和牛顿第二定理，可得：
Figure GDA0004135471660000081
令x1(t)＝y(t)，
Figure GDA0004135471660000082
可得弹簧阻尼系统的系统模型：
Figure GDA0004135471660000083
Figure GDA0004135471660000084
采用零阶保持器法对弹簧阻尼系统的系统模型进行离散化，设定采样时间间隔Ts＝0.1s，则弹簧阻尼系统的系统方程为：
y(k)＝θTφ(k)+e(k) (4)
其中φ(k)为信息向量，并且φ(k)＝[-y(k-1),-y(k-2),u(k-1),u(k-2)]T，θ为参数向量，k表示离散时间，并且θ＝[θ1,θ2,θ3,θ4]T，e(k)为弹簧阻尼系统的不确定噪声，且e(k)有界，即|e(k)|≤γ，γ为大于零的常数。
步骤102，获取弹簧阻尼系统在工作状态下的外加控制力和对应的物块的位移，以确定系统的信息向量。
在预定时间范围内，获取弹簧阻尼系统在工作状态下的外加控制力和对应的物块的位移。
预定时间范围为1至N，N为整数，N的值是预先设置的。
可选的，在弹簧阻尼系统中，利用力传感器测量外加控制力的大小，利用位移传感器测量物块的位移。
将所得的外加控制力和对应的物块的位移的数据代入信息向量的式子φ(k)＝[-y(k-1),-y(k-2),u(k-1),u(k-2)]T中，确定出弹簧阻尼系统的信息向量φ(k)；k的取值范围为1至N，k为整数。
步骤103，根据步骤101确定的弹簧阻尼系统的系统方程和步骤102确定的信息向量，确定参数向量对应的椭球集合和系统故障指示信号的数值。
设置初始化的椭球集合E(0)的中心θc(0)和轴信息矩阵P(0)；设置初始化故障指示信号。
轴信息矩阵是表示椭球集合的形状和大小的一个对称正定矩阵。
将初始化的椭球集合E(0)的中心θc(0)设置为：θc(0)＝[0,0,0,0]T；将初始化的椭球集合E(0)的轴信息矩阵P(0)设置为：P(0)＝δ·In；将初始化的故障指示信号设置为：f(0)＝0。
δ为较大的正数，In为n阶单位矩阵，n为参数向量θ的维数，本实施例中以n＝4为例进行说明。
通过递推获取k时刻(k时刻表示第k个离散时刻，后续简称为k时刻)椭球集合E(k)的中心θc(k)、k时刻椭球集合E(k)的轴信息矩阵P(k)和k时刻弹簧阻尼系统的故障指示信号f(k)。
递推过程如下：
步骤1031，利用k时刻物块的位移y(k)、信息向量φ(k)，k-1时刻的椭球集合的中心θc(k-1)、轴信息矩阵P(k-1)和系统不确定噪声的边界γ，按式(5)和式(6)分别确定k时刻仿射变换后
Figure GDA0004135471660000091
正交的两个平行超平面的
Figure GDA0004135471660000092
坐标，即α(k)和
Figure GDA0004135471660000093
其中
Figure GDA0004135471660000094
为仿射变换后的φ。
Figure GDA0004135471660000095
Figure GDA0004135471660000096
步骤1032，根据α(k)和
Figure GDA0004135471660000097
确定k时刻弹簧阻尼系统的故障指示信号f(k)，具体的：
若α(k)≥1或
Figure GDA0004135471660000098
则故障指示信号f(k)＝1。
若α(k)≤1且
Figure GDA0004135471660000099
则故障指示信号f(k)＝0。
步骤1033，根据α(k)和
Figure GDA00041354716600000910
k时刻物块的位移y(k)、信息向量φ(k)，k-1时刻的椭球集合的中心θc(k-1)、轴信息矩阵P(k-1)和系统不确定噪声的边界γ，更新k时刻椭球集合E(k)的中心θc(k)和轴信息矩阵P(k)。
若α(k)≥1或
Figure GDA00041354716600000911
则按式(7)更新k时刻椭球集合E(k)的中心θc(k)，按式(8)更新k时刻椭球集合E(k)的轴信息矩阵P(k)。
θc(k)＝θc(k-1) (7)
P(k)＝P(k-1) (8)
若α(k)≤1且
Figure GDA00041354716600000912
并且两者同时满足
Figure GDA00041354716600000913
则按式(7)更新k时刻椭球集合E(k)的中心θc(k)，按式(8)更新k时刻椭球集合E(k)的轴信息矩阵P(k)。
n为参数向量θ的维数。
若α(k)≤1且
Figure GDA0004135471660000101
并且两者同时满足
Figure GDA0004135471660000102
则按如下步骤更新k时刻椭球集合E(k)的中心θc(k)和轴信息矩阵P(k)。
S1、根据α(k)和
Figure GDA0004135471660000103
按下式计算得到k时刻两者的平均值μ(k)。
Figure GDA0004135471660000104
S2、若|μ(k)|＞ρ，则按式(10)至式(12)分别计算k时刻仿射变换后
Figure GDA0004135471660000105
沿
Figure GDA0004135471660000106
的中心坐标τ(k)，
Figure GDA0004135471660000107
沿
Figure GDA0004135471660000108
半轴的平方σ(k)，
Figure GDA0004135471660000109
正交于
Figure GDA00041354716600001010
的半轴的平方δ(k)。
Figure GDA00041354716600001011
Figure GDA00041354716600001012
Figure GDA00041354716600001013
其中ρ为大于零的数，可人为进行设定，一般设置为较小的数值，如ρ＝10-6，
Figure GDA00041354716600001014
为中间变量，
Figure GDA00041354716600001015
为仿射变换后的E，
Figure GDA00041354716600001016
为仿射变换后的φ。
若|μ(k)|≤ρ，则按式(13)至式(15)分别计算k时刻仿射变换后
Figure GDA00041354716600001017
沿
Figure GDA00041354716600001018
的中心坐标τ(k)，
Figure GDA00041354716600001019
沿
Figure GDA00041354716600001020
半轴的平方σ(k)，
Figure GDA00041354716600001021
正交
Figure GDA00041354716600001022
的半轴的平方δ(k)。
τ(k)＝0 (13)
σ(k)＝nα2 (14)
Figure GDA00041354716600001023
其中
Figure GDA00041354716600001024
为中间变量。
S3、按式(16)更新k时刻椭球集合E(k)的中心θc(k)，按式(17)更新k时刻椭球集合E(k)的轴信息矩阵P(k)。
Figure GDA00041354716600001025
Figure GDA0004135471660000111
步骤104，根据弹簧阻尼系统参数向量对应的椭球集合，确定椭球集合对应的支持正多胞体集合和支持正多胞体的交集。
根据k时刻弹簧阻尼系统参数向量对应的椭球集合E(k)，按式(18)确定椭球集合对应的支持多胞体集合O(k)：
Figure GDA0004135471660000112
具体地，
Figure GDA0004135471660000113
Figure GDA0004135471660000114
其中，
Figure GDA0004135471660000115
u∈{1,2,…,n}，<·>为内积函数。
根据k时刻椭球集合对应的支持正多胞体集合O(k)，确定支持正多胞体的交集X(k)：
X(k)＝O(1)∩…∩O(k)＝X(k-1)∩O(k) (21)
具体地，按式(22)确定支持正多胞体的交集X(k)。
Figure GDA0004135471660000116
具体地，
Figure GDA0004135471660000117
u∈{1,2,…,n}。
步骤105，根据弹簧阻尼系统的故障指示信号，确定系统的故障状态和故障时间。
若k时刻故障指示信号f(k)＝1，则表明弹簧阻尼系统在k时刻发生故障；
若k时刻故障指示信号f(k)＝0，则表明弹簧阻尼系统在k时刻未发生故障。
步骤106，若系统发生故障，确定弹簧阻尼系统的参数向量θ的具体故障分量。
步骤1061，对k-1时刻支持向量的交集X(k-1)按式(23)进行n-1维扩展，获得测试集合
Figure GDA0004135471660000118
i∈{1,…,n}，即对于第i个测试集合
Figure GDA0004135471660000119
除掉第i维，都进行扩展。
Figure GDA0004135471660000121
其中，
Figure GDA0004135471660000122
为第u维参数θu的最大变化范围，u∈{1,…,n}。
步骤1062，根据k-1时刻的测试集合
Figure GDA0004135471660000123
计算k-1时刻测试集合对应的
Figure GDA0004135471660000124
Figure GDA0004135471660000125
椭球
Figure GDA0004135471660000126
具体地：
根据k-1时刻的测试集合
Figure GDA0004135471660000127
的顶点
Figure GDA0004135471660000128
按如下约束求解
Figure GDA0004135471660000129
椭球
Figure GDA00041354716600001210
的中心
Figure GDA00041354716600001211
i∈{1,…,n}：
Figure GDA00041354716600001212
Qi∈Rn×n为对称正定矩阵，
Figure GDA00041354716600001213
为
Figure GDA00041354716600001214
的第p个顶点。
按式(25)求解
Figure GDA00041354716600001215
椭球
Figure GDA00041354716600001216
的轴信息矩阵
Figure GDA00041354716600001217
Figure GDA00041354716600001218
步骤1063，利用弹簧阻尼系统的系统方程和信息向量，按步骤103中更新参数向量对应的椭球集合的方法更新
Figure GDA00041354716600001219
椭球。
利用k时刻物块的位移y(k)、信息向量φ(k)，k-1时刻的
Figure GDA00041354716600001220
椭球
Figure GDA00041354716600001221
的中心
Figure GDA00041354716600001222
轴信息矩阵Pi t(k-1)和系统不确定噪声的边界γ，按步骤1031和步骤1033的递推更新方法来更新k时刻
Figure GDA00041354716600001232
椭球
Figure GDA00041354716600001223
的中心
Figure GDA00041354716600001224
和轴信息矩阵Pi t(k)，i∈{1,…,n}；利用k+1时刻物块的位移y(k+1)、信息向量φ(k+1)，k时刻的
Figure GDA00041354716600001234
Figure GDA00041354716600001235
椭球
Figure GDA00041354716600001225
的中心
Figure GDA00041354716600001226
轴信息矩阵Pi t(k)和系统不确定噪声的边界γ，按步骤1031和步骤1033的递推更新方法来更新k+1时刻
Figure GDA00041354716600001227
椭球
Figure GDA00041354716600001228
的中心
Figure GDA00041354716600001229
和轴信息矩阵Pi t(k+1)，i∈{1,…,n}；……利用k+L时刻物块的位移y(k+L)、信息向量φ(k+L)，k+L-1时刻的
Figure GDA00041354716600001233
椭球
Figure GDA00041354716600001230
的中心
Figure GDA00041354716600001231
轴信息矩阵Pi t(k+L-1)和系统不确定噪声的边界γ，按步骤1031和步骤1033的递推更新方法来更新k+L时刻
Figure GDA0004135471660000131
椭球
Figure GDA0004135471660000132
的中心
Figure GDA0004135471660000133
和轴信息矩阵Pi t(k+L)，i∈{1,…,n}。
L的值是预先设置的，L≤N-k，L的值根据实际需要确定，比如L＝10。
步骤1064，计算k+L时刻
Figure GDA0004135471660000134
椭球
Figure GDA0004135471660000135
空集指示信号
Figure GDA0004135471660000136
利用k+L时刻物块的位移y(k+L)、信息向量φ(k+L)，k+L-1时刻的
Figure GDA0004135471660000137
椭球
Figure GDA0004135471660000138
的中心
Figure GDA0004135471660000139
轴信息矩阵Pi t(k+L-1)和系统不确定噪声的边界γ，按步骤1031和步骤1032中确定弹簧阻尼系统的故障指示信号的方法来计算k+L时刻
Figure GDA00041354716600001310
Figure GDA00041354716600001311
椭球
Figure GDA00041354716600001312
的空集指示信号
Figure GDA00041354716600001313
步骤1065，根据k+L时刻
Figure GDA00041354716600001314
椭球
Figure GDA00041354716600001315
空集指示信号
Figure GDA00041354716600001316
确定弹簧阻尼系统的参数向量θ的具体故障分量。
具体地：
若对于所有的i，i∈{1,…,n}，当k+L时刻
Figure GDA00041354716600001317
椭球
Figure GDA00041354716600001318
空集指示信号
Figure GDA00041354716600001319
都为1时，则弹簧阻尼系统的参数向量θ中的θi，i∈{1,…,n}，都发生故障；
若对于所有的i，i∈{1,…,n}，当i≠j，j∈{1,…,n}时，k+L时刻
Figure GDA00041354716600001320
椭球
Figure GDA00041354716600001321
空集指示信号
Figure GDA00041354716600001337
都为1，同时当i＝j时，k+L时刻
Figure GDA00041354716600001322
椭球
Figure GDA00041354716600001323
空集指示信号
Figure GDA00041354716600001324
为0，则弹簧阻尼系统的参数向量θ中的θi，i∈{1,…,n}，i≠j都发生故障；
若对于所有的i，i∈{1,…,n}，当i＝j和i＝q，j∈{1,…,n}，q∈{1,…,n}，j≠q时，k+L时刻
Figure GDA00041354716600001325
椭球
Figure GDA00041354716600001326
空集指示信号
Figure GDA00041354716600001327
都为1，同时当i≠j且i≠q时，k+L时刻
Figure GDA00041354716600001328
椭球
Figure GDA00041354716600001329
空集指示信号
Figure GDA00041354716600001330
都为0，则弹簧阻尼系统的参数向量θ中的θi，i∈{1,…,n}，i＝j且i＝q都发生故障；
若对于所有的i，i∈{1,…,n}，当i＝j，j∈{1,…,n}时，k+L时刻
Figure GDA00041354716600001331
椭球
Figure GDA00041354716600001332
空集指示信号
Figure GDA00041354716600001333
为1，同时当i≠j时，k+L时刻
Figure GDA00041354716600001334
椭球
Figure GDA00041354716600001335
空集指示信号
Figure GDA00041354716600001336
都为0，则弹簧阻尼系统的参数向量θ中的θi，i∈{1,…,n}，i＝j发生故障。
步骤107，根据弹簧阻尼系统的参数向量的具体故障分量确定扩展方向，按扩展方向重置支持正多胞体的交集。
根据弹簧阻尼系统的参数向量θ的具体故障分量确定扩展方向，按扩展方向重置k-1时刻支持正多胞体的交集。
若弹簧阻尼系统的参数向量θ中的θi，i∈{1,…,n}发生故障，则按式(26)和式(27)分别更新k-1时刻重置的支持正多胞体的交集Xr(k-1)的第i维参数的上界和下界：
Figure GDA0004135471660000141
Figure GDA0004135471660000142
其中，
Figure GDA0004135471660000143
为第i维参数θi的最大变化范围，i∈{1,…,n}，r为上标，带有此上标的参数和集合表示重置后的参数和集合。
若弹簧阻尼系统的参数向量θ中的θi，i∈{1,…,n}未发生故障，则按式(28)和式(29)分别更新k-1时刻重置的支持正多胞体的交集Xr(k-1)的第i维参数的上界和下界：
Figure GDA0004135471660000144
Figure GDA0004135471660000145
根据k-1时刻重置的支持正多胞体的交集Xr(k-1)每维参数的上界和下界，按式(30)获得k-1时刻按扩展方向重置的支持正多胞体的交集Xr(k-1)：
Figure GDA0004135471660000146
步骤108，根据按扩展方向重置的支持正多胞体的交集，识别弹簧阻尼系统的故障参数。
步骤1081，根据k-1时刻按扩展方向重置的支持正多胞体的交集Xr(k-1)，按步骤1062中计算k-1时刻测试集合
Figure GDA0004135471660000147
对应的
Figure GDA00041354716600001411
椭球
Figure GDA0004135471660000148
的方法计算k-1时刻按扩展方向重置的支持正多胞体的交集Xr(k-1)对应的
Figure GDA0004135471660000149
椭球Er(k-1)。
步骤1082，利用k时刻物块的位移y(k)、信息向量φ(k)，k-1时刻的
Figure GDA00041354716600001410
椭球Er(k-1)的中心θcr(k-1)、轴信息矩阵Pr(k-1)和系统不确定噪声的边界γ，按步骤1031和步骤1033的递推更新方法来重新更新k时刻弹簧阻尼系统参数向量对应的椭球集合E(k)的中心θc(k)和轴信息矩阵P(k)。
步骤1083，利用步骤104重新更新椭球集合对应的支持正多胞体集合和支持正多胞体的交集。
根据k时刻弹簧阻尼系统参数向量对应的椭球集合E(k)，再次更新椭球集合对应的支持多胞体集合O(k)；
根据k时刻椭球集合对应的支持正多胞体集合O(k)，再次更新支持正多胞体的交集X(k)。
步骤1084，根据支持正多胞体的交集X(k)每一维的上界和下界，确定弹簧阻尼系统的故障参数向量。
对于k时刻支持正多胞体的交集X(k)的第i维对应的参数θXi(k)，i∈{1,…,n}，按式(31)计算其参数中心值
Figure GDA0004135471660000151
Figure GDA0004135471660000152
Figure GDA0004135471660000153
为k时刻支持正多胞体的交集X(k)的第i维对应的参数的上界，
Figure GDA0004135471660000154
为k时刻支持正多胞体的交集X(k)的第i维对应的参数的下界，i∈{1,…,n}。
利用k时刻支持正多胞体的交集X(k)的每一维对应的参数的中心值
Figure GDA0004135471660000155
按式(32)确定故障参数向量
Figure GDA0004135471660000156
从而实现k时刻已发生故障的弹簧阻尼系统的故障参数的识别。
Figure GDA0004135471660000157
需要说明的是：本发明实施例提供的面向工业生产过程的弹簧阻尼系统滤波故障诊断方法，诊断弹簧阻尼系统是否发生故障时在弹簧阻尼系统处于工作状态下进行的。
为验证本申请所提出的面向工业生产过程的弹簧阻尼系统滤波故障诊断方法的收敛速度和故障识别时间，进行如下仿真实验：
仿真实验中设置系统噪声为|e(k)|≤0.01，设定的弹簧阻尼系统的相关状态和参数如表1所示，即当时间k到达1001，2001和3001时，我们分别加入故障状态1，故障状态2和故障状态3。
表1弹簧阻尼系统在不同状态下的参数分量的数值
Figure GDA0004135471660000158
基于相同的仿真条件，将本申请提出的一种面向工业生产过程的弹簧阻尼系统滤波故障诊断方法与基于全域放大集员滤波的弹簧阻尼系统故障诊断方法进行对比，得到的故障诊断对比结果如图3至图6所示。
全域放大集员滤波算法可参考“Fault diagnosis based on set membershipidentification using output-error models[J],Int J Adapt Control SignalProcess,2016,30,(2),pp.224–255.”
利用支持正多胞体的交集X(k)的每一维对应的参数的上下界曲线的递推变换过程来展现两种算法对弹簧阻尼系统的故障诊断结果。
以图3为例，在k＝0～4000时刻系统共有4中工作状态：
在k＝0～1000时刻，两种算法所得出的θ1的上下界变化曲线一致，说明在k＝0～1000时刻系统处于正常状态，没有故障出现，两种算法对于正常状态下的系统的参数辨识结果是一致的。
在k＝1001～2000时刻，此时系统有故障发生，系统参数发生变化，但是θ1并没有发生变化，因此本申请所提出的按扩展方向进行重置的弹簧阻尼系统滤波故障诊断方法可以基于故障分量，得出θ1并没有发生变化，无须在θ1的方向上对支持正多胞体的交集进行重置，因此θ1的上下界变化曲线继续保持上一状态(即k＝1～1000时刻)的收缩趋势，上下界之间范围越小，就可以得到更精确的系统的参数值；而全域放大的算法，在系统的每一个状态都对支持正多胞体的交集在每一个参数分量的方向上都进行了重置，因此在k＝1001～2000时刻，尽管θ1没有发生故障，但是依然对支持正多胞体的交集在θ1方向上进行了重置，直接导致的结果是需要再次对θ1的上下界从一个很大的范围进行收缩，因此该算法下的θ1的上下界范围比按扩展方向进行重置的弹簧阻尼系统滤波故障诊断方法大，进而算法的收敛速度更慢，识别出系统的故障参数所需时间越长。
在k＝2001～3000时刻，系统发生故障，并且参数分量θ1发生故障，本专利所提出的按扩展方向进行重置的弹簧阻尼系统滤波故障诊断方法基于故障分量，得出θ1发生故障，在θ1的方向上对支持正多胞体的交集进行了重置，同时全域放大的故障诊断方法也对支持正多胞体的交集在θ1的方向上进行了重置，但仍然可以看出，本专利所提出的按扩展方向进行重置的弹簧阻尼系统滤波故障诊断方法下的θ1的上下界之间的范围收缩得更快。
在k＝3001～4000时刻，系统发生故障，但是θ1并没有发生变化，此时两种算法的故障诊断效果与k＝1001～2000时刻一致，故障诊断分析也一致。
因此，根据图3至图6，可以得出以下结论：
(1)基于两种算法，X(k)的对应的参数的上下界在整个故障诊断过程中被重置了三次，说明此弹簧阻尼系统在整个工作过程中发生了三次故障。
(2)在第一个状态下两种算法得出的X(k)有着相同的参数的上下界，并且通过表2可得，两种算法有着相同的参数辨识结果。说明第一个状态为正常状态，并且两种算法都有着较高的参数辨识精度。
(3)精确分析可得，参数分量的上下界在k＝1001，2001和3001时刻分别进行了重置，说明两种算法都能快速及时地检测到系统的故障，确定出系统故障时间。
(4)面向工业生产过程的弹簧阻尼系统滤波故障诊断方法在故障状态1中重置了第三维参数和第四维参数，在故障状态2中重置了第一维参数和第三维参数，在故障状态3中重置了第二维参数、第三维参数和第四维参数。因此，可以观察得出第三个和第四个参数分量在故障状态1中被进行了故障隔离，第一个和第三个参数分量在故障状态2中被进行了故障隔离，第二个，第三个和第四个参数分量在故障状态4中被进行了故障隔离，这些具体故障分量的诊断结果与系统的真实故障状态相一致，说明该算法能够进行快速的故障隔离。
(5)根据故障隔离的结果，面向工业生产过程的弹簧阻尼系统滤波故障诊断方法能够根据具体故障分量定向重置正多胞体的交集，相比于基于全域放大集员滤波的弹簧阻尼系统故障诊断方法，该方法有着更小的上下界区间，相应地，该方法对于弹簧阻尼系统的故障收敛速度更快，能够更快地识别出系统的故障参数。
(6)表2中展示了两种算法的最终故障辨识结果，说明两种算法都能够精确地识别出弹簧阻尼系统的故障。
表2两种算法的故障识别终值
Figure GDA0004135471660000171
需要说明的是，表2中
Figure GDA0004135471660000172
表示向工业生产过程的弹簧阻尼系统滤波故障诊断方法所得出的参数识别结果，
Figure GDA0004135471660000173
表示基于全域放大集员滤波的弹簧阻尼系统故障诊断方法所得出的参数识别结果。对比表1可知，这两种算法在识别精度上都比较高。但是结合图3-图6可知，本申请所提出的工业生产过程的弹簧阻尼系统滤波故障诊断方法能更快地识别出故障。
综上所述，本申请针对工业生产过程中弹簧阻尼系统故障诊断问题，提出了一类面向工业生产过程的弹簧阻尼系统滤波故障诊断方法，通过建立弹簧阻尼系统的离散系统模型，确定弹簧阻尼系统的系统方程；获取弹簧阻尼系统在工作状态下的外加控制力和物块的位移，确定系统的信息向量；根据弹簧阻尼系统的系统方程和信息向量，确定参数向量对应的椭球集合和系统故障指示信号的数值；根据弹簧阻尼系统参数向量对应的椭球集合，确定椭球集合对应的支持正多胞体集合和支持正多胞体的交集；根据弹簧阻尼系统的故障指示信号，确定系统的故障状态和故障时间；若系统发生故障，确定弹簧阻尼系统的参数向量的具体故障分量；根据弹簧阻尼系统的参数向量的具体故障分量，按扩展方向重置支持正多胞体的交集；根据按扩展方向重置的支持正多胞体的交集，识别弹簧阻尼系统的故障参数；解决了系统噪声不确定下的弹簧阻尼系统的快速准确实时的故障诊断，同时相比于基于全域放大集员滤波的故障诊断方法，面向工业生产过程的弹簧阻尼系统滤波故障诊断方法在弹簧阻尼系统的故障诊断中有着收敛速度快，故障识别时间短的优点。
本发明实施例中的部分步骤，可以利用软件实现，相应的软件程序可以存储在可读取的存储介质中，如光盘或硬盘等。
以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种弹簧阻尼系统滤波故障诊断方法，其特征在于，所述方法包括：
根据弹簧阻尼系统的系统方程和信息向量确定弹簧阻尼系统的参数向量对应的椭球集合，进而确定椭球集合对应的支持正多胞体集合和支持正多胞体的交集；
在系统发生故障时，根据测试集合及其
Figure FDA0004135471650000011
-John椭球确定弹簧阻尼系统的参数向量的具体故障分量，继而根据弹簧阻尼系统的参数向量的具体故障分量确定扩展方向，按扩展方向重置支持正多胞体的交集；
根据按扩展方向重置的支持正多胞体的交集，识别弹簧阻尼系统的参数向量中的故障参数；
所述方法包括：
步骤101，建立弹簧阻尼系统的离散系统模型，确定弹簧阻尼系统的系统方程；
y(k)＝θTφ(k)+e(k)
其中，y(k)为k时刻系统的输出，θ＝[θ1,θ2,θ3,θ4]T为参数向量，k表示离散时间，φ(k)为信息向量，e(k)为弹簧阻尼系统的不确定噪声，且e(k)有界，即|e(k)|≤γ，γ为大于零的常数；
步骤102，获取弹簧阻尼系统在工作状态下的外加控制力和对应的物块的位移，以确定系统的信息向量φ(k)；
步骤103，根据步骤101确定的弹簧阻尼系统的系统方程和步骤102确定的信息向量φ(k)，确定参数向量θ对应的椭球集合和系统故障指示信号f(k)的数值；
步骤104，根据步骤103确定的弹簧阻尼系统参数向量θ对应的椭球集合，确定椭球集合对应的支持正多胞体集合O(k)和支持正多胞体的交集；
步骤105，根据步骤103确定的弹簧阻尼系统的故障指示信号f(k)，确定系统的故障状态和故障时间；
步骤106，若系统发生故障，确定弹簧阻尼系统的参数向量θ的具体故障分量；
步骤107，根据步骤106确定的弹簧阻尼系统的参数向量的具体故障分量确定扩展方向，按扩展方向重置支持正多胞体的交集；
步骤108，根据步骤107中按扩展方向重置的支持正多胞体的交集，识别弹簧阻尼系统的故障参数；
所述步骤104，根据步骤103确定的弹簧阻尼系统参数向量θ对应的椭球集合，确定椭球集合对应的支持正多胞体集合O(k)和支持正多胞体的交集，包括：
Figure FDA0004135471650000021
Figure FDA0004135471650000022
Figure FDA0004135471650000023
其中，
Figure FDA0004135471650000024
表示k时刻支持正多胞体集合O(k)的参数的上界，
Figure FDA0004135471650000025
表示k时刻支持正多胞体集合O(k)的参数的下界，
Figure FDA0004135471650000026
u∈{1,2,…,n}，<·>为内积函数，n为参数向量θ的维数，k时刻表示第k个离散时刻；
根据支持正多胞体集合O(k)，确定k时刻支持正多胞体的交集：
X(k)＝O(1)∩…∩O(k)＝X(k-1)∩O(k)
即
Figure FDA0004135471650000027
Figure FDA0004135471650000028
Figure FDA0004135471650000029
其中，
Figure FDA00041354716500000210
表示k时刻支持正多胞体的交集X(k)的参数的上界，
Figure FDA00041354716500000211
表示k时刻支持正多胞体的交集X(k)的参数的下界。
2.根据权利要求1所述的方法，其特征在于，所述步骤106，若系统发生故障，确定弹簧阻尼系统的参数向量θ的具体故障分量，包括：
对k-1时刻支持正多胞体的交集X(k-1)进行n-1维扩展，获得测试集合
Figure FDA00041354716500000223
i∈{1,…,n}，即对于第i个测试集合
Figure FDA00041354716500000224
除掉第i维，都进行扩展；
根据所述k-1时刻的测试集合
Figure FDA00041354716500000225
计算k-1时刻测试集合对应的
Figure FDA00041354716500000212
-John椭球
Figure FDA00041354716500000226
利用弹簧阻尼系统的系统方程和信息向量，更新k-1时刻的
Figure FDA00041354716500000213
-John椭球
Figure FDA00041354716500000214
得到k时刻的
Figure FDA00041354716500000215
-John椭球
Figure FDA00041354716500000216
计算k+L时刻
Figure FDA00041354716500000217
-John椭球
Figure FDA00041354716500000218
空集指示信号
Figure FDA00041354716500000219
根据所述k+L时刻
Figure FDA00041354716500000220
-John椭球
Figure FDA00041354716500000221
空集指示信号
Figure FDA00041354716500000222
确定弹簧阻尼系统的参数向量θ的具体故障分量。
3.根据权利要求2所述的方法，其特征在于，所述步骤107，根据步骤106确定的弹簧阻尼系统的参数向量的具体故障分量确定扩展方向，按扩展方向重置支持正多胞体的交集，包括：
若弹簧阻尼系统的参数向量θ中的θi发生故障，按下式更新k-1时刻重置的支持正多胞体的交集Xr(k-1)的第i维参数的上界：
Figure FDA0004135471650000031
若弹簧阻尼系统的参数向量θ中的θi发生故障，按下式更新k-1时刻重置的支持正多胞体的交集Xr(k-1)的第i维参数的下界：
Figure FDA0004135471650000032
若弹簧阻尼系统的参数向量θ中的θi未发生故障，按下式更新k-1时刻重置的支持正多胞体的交集Xr(k-1)的第i维参数的上界：
Figure FDA0004135471650000033
若弹簧阻尼系统的参数向量θ中的θi未发生故障，按下式更新k-1时刻重置的支持正多胞体的交集Xr(k-1)的第i维参数的下界：
Figure FDA0004135471650000034
根据所述k-1时刻支持正多胞体的交集Xr(k-1)每维参数的上界和下界，按下式获得k-1时刻按扩展方向重置的支持正多胞体的交集Xr(k-1)：
Figure FDA0004135471650000035
其中，
Figure FDA0004135471650000036
为第i维参数θi的最大变化范围，i∈{1,…,n}，r为上标，带有此上标的参数和集合表示重置后的参数和集合。
4.根据权利要求3所述的方法，其特征在于，所述根据步骤107中按扩展方向重置的支持正多胞体的交集，识别弹簧阻尼系统的故障参数，包括：
根据所述k-1时刻按扩展方向重置的支持正多胞体的交集Xr(k-1)，计算k-1时刻按扩展方向重置的支持正多胞体的交集Xr(k-1)对应的
Figure FDA0004135471650000037
-John椭球Er(k-1)；
利用k时刻物块的位移y(k)、信息向量φ(k)，k-1时刻的
Figure FDA0004135471650000038
-John椭球Er(k-1)的中心θcr(k-1)、轴信息矩阵Pr(k-1)和系统不确定噪声的边界γ，重新更新k时刻弹簧阻尼系统参数向量对应的椭球集合E(k)的中心θc(k)和轴信息矩阵P(k)；
根据所述重新更新的椭球集合E(k)，重新更新k时刻椭球集合E(k)对应的支持正多胞体集合O(k)和支持正多胞体的交集X(k)；
根据所述支持正多胞体的交集X(k)每一维的上界和下界，确定弹簧阻尼系统的故障参数向量
Figure FDA00041354716500000413
Figure FDA0004135471650000041
其中
Figure FDA0004135471650000042
为k时刻支持正多胞体的交集X(k)的第i维对应的参数的上界，
Figure FDA0004135471650000043
为k时刻支持正多胞体的交集X(k)的第i维对应的参数的下界。
5.根据权利要求4所述的方法，其特征在于，所述步骤103，根据步骤101确定的弹簧阻尼系统的系统方程和步骤102确定的信息向量φ(k)，确定参数向量θ对应的椭球集合和系统故障指示信号f(k)的数值，包括：
按下面两式确定k时刻仿射变换后
Figure FDA0004135471650000044
正交的两个平行超平面的
Figure FDA0004135471650000045
坐标，即α(k)和
Figure FDA0004135471650000046
其中
Figure FDA0004135471650000047
为仿射变换后的φ：
Figure FDA0004135471650000048
Figure FDA0004135471650000049
若α(k)≥1或
Figure FDA00041354716500000410
则故障指示信号f(k)＝1，表明弹簧阻尼系统在k时刻发生故障；
θc(k)＝θc(k-1)，
P(k)＝P(k-1)；
若α(k)≤1且
Figure FDA00041354716500000411
则故障指示信号f(k)＝0，表明弹簧阻尼系统在k时刻未发生故障。
6.根据权利要求5所述的方法，其特征在于，若
Figure FDA00041354716500000412
按下面两式更新k时刻椭球集合E(k)的中心θc(k)和轴信息矩阵P(k)：
θc(k)＝θc(k-1)，
P(k)＝P(k-1)；
若两者还同时满足
Figure FDA0004135471650000051
则
Figure FDA0004135471650000052
在上述条件下，若|μ(k)|＞ρ，则
Figure FDA0004135471650000053
Figure FDA0004135471650000054
Figure FDA0004135471650000055
Figure FDA0004135471650000056
Figure FDA0004135471650000057
在上述条件下，若|μ(k)|≤ρ，则
Figure FDA0004135471650000058
τ(k)＝0，
σ(k)＝nα2，
Figure FDA0004135471650000059
之后按下面两式更新k时刻椭球集合E(k)的中心θc(k)和轴信息矩阵P(k)
Figure FDA00041354716500000510
Figure FDA00041354716500000511
其中n为参数向量θ的维数，μ(k)为α(k)和
Figure FDA00041354716500000512
的平均值，ρ为大于零的数，设定ρ＝10-6，τ(k)为k时刻仿射变换后
Figure FDA00041354716500000513
沿
Figure FDA00041354716500000514
的中心坐标，σ(k)为
Figure FDA00041354716500000515
沿
Figure FDA00041354716500000516
半轴的平方，δ(k)为
Figure FDA00041354716500000517
正交于
Figure FDA00041354716500000518
的半轴的平方，ε(k)、b(k)、α为中间变量，
Figure FDA00041354716500000519
为仿射变换后的E，
Figure FDA00041354716500000520
为仿射变换后的φ。
7.根据权利要求6所述的方法，其特征在于，所述步骤102，获取弹簧阻尼系统在工作状态下的外加控制力和对应的物块的位移，以确定系统的信息向量φ(k)，包括：
在预定时间范围内，获取弹簧阻尼系统在工作状态下的外加控制力和对应的物块的位移；
将所得的外加控制力和对应的物块的位移的数据代入下式中：
φ(k)＝[-y(k-1),-y(k-2),u(k-1),u(k-2)]T
确定出弹簧阻尼系统的信息向量φ(k)；k的取值范围为1至N，k为整数。
8.一种弹簧阻尼系统滤波故障诊断系统，其特征在于，所述系统采用权利要求1-7任一所述的方法对弹簧阻尼系统进行故障诊断。
CN202010272258.2A 2020-04-09 2020-04-09 一种面向工业生产过程的弹簧阻尼系统滤波故障诊断方法 Active CN111597647B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010272258.2A CN111597647B (zh) 2020-04-09 2020-04-09 一种面向工业生产过程的弹簧阻尼系统滤波故障诊断方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010272258.2A CN111597647B (zh) 2020-04-09 2020-04-09 一种面向工业生产过程的弹簧阻尼系统滤波故障诊断方法

Publications (2)

Publication Number Publication Date
CN111597647A CN111597647A (zh) 2020-08-28
CN111597647B true CN111597647B (zh) 2023-04-25

Family

ID=72181898

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010272258.2A Active CN111597647B (zh) 2020-04-09 2020-04-09 一种面向工业生产过程的弹簧阻尼系统滤波故障诊断方法

Country Status (1)

Country Link
CN (1) CN111597647B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112149886B (zh) * 2020-09-07 2023-10-31 江南大学 一种基于多维空间滤波的四容水箱系统状态估计方法
CN112305418B (zh) * 2020-10-13 2021-09-28 江南大学 一种基于混合噪声双重滤波的电机系统故障诊断方法
CN112883508B (zh) * 2021-01-22 2024-03-08 江南大学 一种基于平行空间滤波的弹簧阻尼系统状态估计方法
CN113359667B (zh) * 2021-06-04 2022-07-22 江南大学 一种基于凸空间滤波的工业系统故障诊断方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
DE102004025574A1 (de) * 2004-05-25 2006-01-19 Abb Research Ltd. Verfahren zur Fehlererkennung in einem industriellen Prozess
WO2017128455A1 (zh) * 2016-01-25 2017-08-03 合肥工业大学 一种基于广义多核支持向量机的模拟电路故障诊断方法
CN109932598A (zh) * 2019-03-25 2019-06-25 江南大学 一种不确定噪声扰动下Buck变换器的故障检测方法
CN110259647A (zh) * 2019-06-21 2019-09-20 江南大学 一种基于正多胞体滤波的风电机组缓变型故障诊断方法
CN110489891A (zh) * 2019-08-23 2019-11-22 江南大学 一种基于多胞空间滤波的工业过程时变参数估计方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
DE102004025574A1 (de) * 2004-05-25 2006-01-19 Abb Research Ltd. Verfahren zur Fehlererkennung in einem industriellen Prozess
WO2017128455A1 (zh) * 2016-01-25 2017-08-03 合肥工业大学 一种基于广义多核支持向量机的模拟电路故障诊断方法
CN109932598A (zh) * 2019-03-25 2019-06-25 江南大学 一种不确定噪声扰动下Buck变换器的故障检测方法
CN110259647A (zh) * 2019-06-21 2019-09-20 江南大学 一种基于正多胞体滤波的风电机组缓变型故障诊断方法
CN110489891A (zh) * 2019-08-23 2019-11-22 江南大学 一种基于多胞空间滤波的工业过程时变参数估计方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
V. Reppa et al.Fault diagnosis based on set membership identification using output-error models.《INTERNATIONAL JOURNAL OF ADAPTIVE CONTROL AND SIGNAL PROCESSING》.2015,第30卷(第2期),全文. *

Also Published As

Publication number Publication date
CN111597647A (zh) 2020-08-28

Similar Documents

Publication Publication Date Title
CN111597647B (zh) 一种面向工业生产过程的弹簧阻尼系统滤波故障诊断方法
Huang et al. Sensor fault diagnosis for structural health monitoring based on statistical hypothesis test and missing variable approach
CN103197663B (zh) 一种故障预测方法及系统
Kourehli et al. Structural damage detection using incomplete modal data and incomplete static response
CN112001110B (zh) 一种基于振动信号空间时时递归图卷积神经网络的结构损伤识别监测方法
CN111639678B (zh) 一种基于集成神经网络的ins/gps组合导航故障检测与诊断方法
CN112633339A (zh) 轴承故障智能诊断方法、诊断系统、计算机设备及介质
CN104298870A (zh) 一种移动荷载下简支梁损伤和移动力同时识别方法
US8566057B2 (en) Method for self-adjustment of a triaxial acceleration sensor during operation, and sensor system having a three-dimensional acceleration sensor
Liu et al. A novel strategy for response and force reconstruction under impact excitation
CN115455353A (zh) 一种基于非线性时域滤波的在线参数辨识方法
Lakshmi et al. Structural damage detection using ARMAX time series models and cepstral distances
CN111859250A (zh) 一种环境及受迫激励下快速模态参数准确性评估方法、系统及存储介质
Liu et al. Performance degradation assessment for coaxial bearings using kernel JADE and two-class model
CN115070765B (zh) 一种基于变分推断的机器人状态估计方法及系统
CN115218927B (zh) 基于二次卡尔曼滤波的无人机imu传感器故障检测方法
CN113392524B (zh) 一种传感器的漂移诊断方法、装置、电子设备及存储介质
CN114139461B (zh) 一种航天器系统不确定性预测方法、装置及存储介质
Feizi et al. Identifying damage location under statistical pattern recognition by new feature extraction and feature analysis methods
CN114234970A (zh) 一种运动载体姿态实时测量方法、装置
Döhler et al. Efficient computation of minmax tests for fault isolation and their application to structural damage localization
CN110738248A (zh) 状态感知数据特征提取方法及装置、系统性能评估方法
US20130179132A1 (en) Analysis Method, Apparatus and Software for a System With Frequency Dependent Materials
CN117932199B (zh) 一种时域声源定位方法
CN117807373A (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