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

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

Info

Publication number
CN111597647A
CN111597647A CN202010272258.2A CN202010272258A CN111597647A CN 111597647 A CN111597647 A CN 111597647A CN 202010272258 A CN202010272258 A CN 202010272258A CN 111597647 A CN111597647 A CN 111597647A
Authority
CN
China
Prior art keywords
spring damping
damping system
time
fault
intersection
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.)
Granted
Application number
CN202010272258.2A
Other languages
English (en)
Other versions
CN111597647B (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 RE-DDA0002581961280000011
椭球确定弹簧阻尼系统的参数向量θ的具体故障分量，继而根据具体故障分量得到扩展方向，按扩展方向重置支持正多胞体的交集，相对于现有故障检测方法中的全域扩展，本申请所提出的面向工业生产过程的弹簧阻尼系统滤波故障诊断方法收敛速度更快、可以更快的识别出故障参数、从而实时性能更好。

Description

一种面向工业生产过程的弹簧阻尼系统滤波故障诊断方法
技术领域
本发明涉及一种面向工业生产过程的弹簧阻尼系统滤波故障诊断方法，属于弹簧阻尼系统故障诊断技术领域。
背景技术
弹簧阻尼系统是工业生产过程涉及的机械振动系统当中普遍使用的一类结构，在工业生产过程中主要起到吸收和耗散生产过程能量的作用，其吸收和耗散能量的大小，关系到生产过程的安全稳定，所以其系统的安全性将会对整个系统的稳定运行产生重大的影响。一旦弹簧阻尼系统发生故障而未被及时诊断，那么将会影响整个系统的运行。因此，对弹簧阻尼系统进行故障诊断的研究对工业生产过程有着重要的现实意义。
常见的故障诊断方法主要分为基于模型的故障诊断方法和基于数据的故障诊断方法两大类。其中基于数据的故障诊断方法需以大量的数据分析为前提，很难做到实时故障诊断。基于模型的故障诊断方法以系统解析模型为基础，能够实现实时在线的故障诊断；但是一般的基于模型的故障诊断方法中要求系统噪声符合一定概率分布，然而实际应用中，弹簧阻尼系统常受各种外部因素的影响，噪声干扰因素不确定，导致其无法满足这一要求。
集员滤波方法只要求系统噪声有界，对噪声的概率分布没有限定，因此基于集员估计的故障诊断方法能够有效处理噪声不确定系统的故障诊断问题。集员估计方法中是基于可行集的状态来确定系统的故障状态。若检测到可行集为空，则认为系统存在故障；若系统可行集不为空，则认为系统没有故障。集员估计方法中根据近似可行集的空间形状和包围方法的不同，可将其划分为椭球体、全对称多胞体、正多胞体等等。
现有的基于集员滤波理论的故障诊断方法中，在检测出系统发生故障后，若要通过重置参数集合的方法来进一步识别故障，一般是在对参数集合进行全域扩展下进行的，然而全域扩展下的参数集合在参数向量的各个维度都需要再次进行收缩，所以存在着收敛速度慢、故障识别时间长，实时性较差的问题。
发明内容
为了解决目前存在的问题，本发明提供了一种面向工业生产过程的弹簧阻尼系统滤波故障诊断方法，所述方法包括：
根据弹簧阻尼系统的系统方程和信息向量确定弹簧阻尼系统的参数向量对应的椭球集合，进而确定椭球集合对应的支持正多胞体集合和支持正多胞体的交集；
在系统发生故障时，根据测试集合及其
Figure BDA0002443506160000022
-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 BDA0002443506160000021
Figure BDA0002443506160000031
Figure BDA0002443506160000032
其中，
Figure BDA0002443506160000033
表示k时刻支持正多胞体集合O(k)的参数的上界，
Figure BDA0002443506160000034
表示k时刻支持正多胞体集合O(k)的参数的下界，
Figure BDA00024435061600000327
<·>为内积函数，n为参数向量θ的维数，k时刻表示第k个离散时刻；
根据支持正多胞体集合O(k)，确定k时刻支持正多胞体的交集：
X(k)＝O(1)∩…∩O(k)＝X(k-1)∩O(k)
即
Figure BDA0002443506160000036
Figure BDA0002443506160000037
Figure BDA0002443506160000038
其中，
Figure BDA0002443506160000039
表示k时刻支持正多胞体的交集X(k)的参数的上界，
Figure BDA00024435061600000310
表示k时刻支持正多胞体的交集X(k)的参数的下界。
可选的，所述步骤106，若系统发生故障，确定弹簧阻尼系统的参数向量θ的具体故障分量，包括：
对k-1时刻支持正多胞体的交集X(k-1)进行n-1维扩展，获得测试集合
Figure BDA00024435061600000311
Figure BDA00024435061600000312
即对于第i个测试集合
Figure BDA00024435061600000325
除掉第i维，都进行扩展；
根据所述k-1时刻的测试集合
Figure BDA00024435061600000326
计算k-1时刻测试集合对应的
Figure BDA00024435061600000313
-John椭球
Figure BDA00024435061600000314
利用弹簧阻尼系统的系统方程和信息向量，更新k-1时刻的
Figure BDA00024435061600000324
-John椭球
Figure BDA00024435061600000315
得到k时刻的
Figure BDA00024435061600000316
John椭球
Figure BDA00024435061600000317
计算所述k+L时刻
Figure BDA00024435061600000318
-John椭球
Figure BDA00024435061600000319
空集指示信号
Figure BDA00024435061600000320
根据所述k+L时刻
Figure BDA00024435061600000321
-John椭球
Figure BDA00024435061600000322
空集指示信号
Figure BDA00024435061600000323
确定弹簧阻尼系统的参数向量θ的具体故障分量。
可选的，所述步骤107，根据步骤106确定的弹簧阻尼系统的参数向量的具体故障分量确定扩展方向，按扩展方向重置支持正多胞体的交集，包括：
若弹簧阻尼系统的参数向量θ中的θi发生故障，按下式更新k-1时刻重置的支持正多胞体的交集Xr(k-1)的第i维参数的上界：
Figure BDA0002443506160000041
若弹簧阻尼系统的参数向量θ中的θi发生故障，按下式更新k-1时刻重置的支持正多胞体的交集Xr(k-1)的第i维参数的下界：
Figure BDA0002443506160000042
若弹簧阻尼系统的参数向量θ中的θi未发生故障，按下式更新k-1时刻重置的支持正多胞体的交集Xr(k-1)的第i维参数的上界：
Figure BDA0002443506160000043
若弹簧阻尼系统的参数向量θ中的θi未发生故障，按下式更新k-1时刻重置的支持正多胞体的交集Xr(k-1)的第i维参数的下界：
Figure BDA0002443506160000044
根据所述k-1时刻支持正多胞体的交集Xr(k-1)每维参数的上界和下界，按下式获得k-1时刻按扩展方向重置的支持正多胞体的交集Xr(k-1)：
Figure BDA0002443506160000045
其中，
Figure BDA0002443506160000046
为第i维参数θi的最大变化范围，i∈{1,…,n}，r为上标，带有此上标的参数和集合表示重置后的参数和集合。
可选的，所述根据步骤107中按扩展方向重置的支持正多胞体的交集，识别弹簧阻尼系统的故障参数，包括：
根据所述k-1时刻按扩展方向重置的支持正多胞体的交集Xr(k-1)，计算k-1时刻按扩展方向重置的支持正多胞体的交集Xr(k-1)对应的
Figure BDA0002443506160000047
-John椭球Er(k-1)；
利用k时刻物块的位移y(k)、信息向量φ(k)，k-1时刻的
Figure BDA0002443506160000048
-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 BDA0002443506160000051
Figure BDA0002443506160000052
其中
Figure BDA0002443506160000053
Figure BDA0002443506160000054
为k时刻支持正多胞体的交集X(k)的第i维对应的参数的上界，
Figure BDA0002443506160000055
为k时刻支持正多胞体的交集X(k)的第i维对应的参数的下界。
可选的，所述步骤103，根据步骤101确定的弹簧阻尼系统的系统方程和步骤102确定的信息向量φ(k)，确定参数向量θ对应的椭球集合和系统故障指示信号f(k)的数值，包括：
按下面两式确定k时刻仿射变换后
Figure BDA0002443506160000056
正交的两个平行超平面的
Figure BDA0002443506160000057
坐标，即α(k)和
Figure BDA0002443506160000058
其中
Figure BDA0002443506160000059
为仿射变换后的φ：
Figure BDA00024435061600000510
Figure BDA00024435061600000511
若α(k)≥1或
Figure BDA00024435061600000512
则故障指示信号f(k)＝1，表明弹簧阻尼系统在k时刻发生故障；
θc(k)＝θc(k-1)，
P(k)＝P(k-1)；
若α(k)≤1且
Figure BDA00024435061600000513
则故障指示信号f(k)＝0，表明弹簧阻尼系统在k时刻未发生故障。
可选的，若
Figure BDA00024435061600000514
按下面两式更新k时刻椭球集合E(k)的中心θc(k)和轴信息矩阵P(k)：
θc(k)＝θc(k-1)，
P(k)＝P(k-1)；
若两者还同时满足
Figure BDA00024435061600000515
则
Figure BDA00024435061600000516
在上述条件下，若|μ(k)|＞ρ，则
Figure BDA0002443506160000061
Figure BDA0002443506160000062
Figure BDA0002443506160000063
Figure BDA0002443506160000064
Figure BDA0002443506160000065
在上述条件下，若|μ(k)|≤ρ，则
Figure BDA0002443506160000066
τ(k)＝0，
σ(k)＝nα2，
Figure BDA0002443506160000067
之后按下面两式更新k时刻椭球集合E(k)的中心θc(k)和轴信息矩阵P(k)
Figure BDA0002443506160000068
Figure BDA0002443506160000069
其中n为参数向量θ的维数，μ(k)为α(k)和
Figure BDA00024435061600000610
的平均值，ρ为大于零的数，设定ρ＝10-6，τ(k)为k时刻仿射变换后
Figure BDA00024435061600000611
沿
Figure BDA00024435061600000612
的中心坐标，σ(k)为
Figure BDA00024435061600000613
沿
Figure BDA00024435061600000614
半轴的平方，δ(k)为
Figure BDA00024435061600000615
正交于
Figure BDA00024435061600000616
的半轴的平方，ε(k)、b(k)、α为中间变量，
Figure BDA00024435061600000617
为仿射变换后的E，
Figure BDA00024435061600000618
为仿射变换后的φ。
可选的，所述步骤102，获取弹簧阻尼系统在工作状态下的外加控制力和对应的物块的位移，以确定系统的信息向量φ(k)，包括：
在预定时间范围内，获取弹簧阻尼系统在工作状态下的外加控制力和对应的物块的位移；
将所得的外加控制力和对应的物块的位移的数据代入下式中：
φ(k)＝[-y(k-1),-y(k-2),u(k-1),u(k-2)]T
确定出弹簧阻尼系统的信息向量φ(k)；k的取值范围为1至N，k为整数。
本申请的另一个方面还提供一种弹簧阻尼系统滤波故障诊断系统，所述系统采用上述面向工业生产过程的弹簧阻尼系统滤波故障诊断方法对弹簧阻尼系统进行故障诊断。
本发明有益效果是：
本申请通过弹簧阻尼系统故障指示信号的数值，确定系统的故障状态，在检测出系统发生故障时，根据测试集合及其-John椭球确定弹簧阻尼系统的参数向量θ的具体故障分量，继而根据具体故障分量得到扩展方向，按扩展方向重置支持正多胞体的交集，相对于现有故障检测方法中的全域扩展，本申请所提出的面向工业生产过程的弹簧阻尼系统滤波故障诊断方法收敛速度更快、可以更快的识别出故障参数、从而实时性能更好。
附图说明
为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1是本发明一个实施例中公开的一种面向工业生产过程的弹簧阻尼系统滤波故障诊断方法流程图。
图2是一种弹簧系统结构图。
图3是本发明一个实施例中公开的两种算法下的一种弹簧系统的参数分量θ1的上下界变化曲线图。
图4是本发明一个实施例中公开的两种算法下的一种弹簧系统的参数分量θ2的上下界变化曲线图。
图5是本发明一个实施例中公开的两种算法下的一种弹簧系统的参数分量θ3的上下界变化曲线图。
图6是本发明一个实施例中公开的两种算法下的一种弹簧系统的参数分量θ4的上下界变化曲线图。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地详细描述。
实施例一：
本实施例提供一种面向工业生产过程的弹簧阻尼系统滤波故障诊断方法，参见图1，所述方法包括：
步骤101，建立弹簧阻尼系统的离散系统模型，确定弹簧阻尼系统的系统方程。
弹簧阻尼系统的结构图如图2所示，其中M为物块的质量，C为摩擦系数，K是弹簧的弹性系数，u(t)为t时刻的外加控制力，是系统的输入，y(t)为t时刻的物块的位移，是系统的输出。
根据此模型结构和牛顿第二定理，可得：
Figure BDA0002443506160000081
令x1(t)＝y(t)，
Figure BDA0002443506160000082
可得弹簧阻尼系统的系统模型：
Figure BDA0002443506160000083
Figure BDA0002443506160000084
采用零阶保持器法对弹簧阻尼系统的系统模型进行离散化，设定采样时间间隔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 BDA0002443506160000091
正交的两个平行超平面的
Figure BDA0002443506160000092
坐标，即α(k)和
Figure BDA0002443506160000093
其中
Figure BDA0002443506160000094
为仿射变换后的φ。
Figure BDA0002443506160000095
Figure BDA0002443506160000096
步骤1032，根据α(k)和
Figure BDA0002443506160000097
确定k时刻弹簧阻尼系统的故障指示信号f(k)，具体的：
若α(k)≥1或
Figure BDA0002443506160000098
则故障指示信号f(k)＝1。
若α(k)≤1且
Figure BDA0002443506160000099
则故障指示信号f(k)＝0。
步骤1033，根据α(k)和
Figure BDA00024435061600000910
k时刻物块的位移y(k)、信息向量φ(k)，k-1时刻的椭球集合的中心θc(k-1)、轴信息矩阵P(k-1)和系统不确定噪声的边界γ，更新k时刻椭球集合E(k)的中心θc(k)和轴信息矩阵P(k)。
若α(k)≥1或
Figure BDA00024435061600000911
则按式(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 BDA0002443506160000101
并且两者同时满足
Figure BDA0002443506160000102
则按式(7)更新k时刻椭球集合E(k)的中心θc(k)，按式(8)更新k时刻椭球集合E(k)的轴信息矩阵P(k)。
n为参数向量θ的维数。
若α(k)≤1且
Figure BDA0002443506160000103
并且两者同时满足
Figure BDA0002443506160000104
则按如下步骤更新k时刻椭球集合E(k)的中心θc(k)和轴信息矩阵P(k)。
S1、根据α(k)和
Figure BDA0002443506160000105
按下式计算得到k时刻两者的平均值μ(k)。
Figure BDA0002443506160000106
S2、若|μ(k)|＞ρ，则按式(10)至式(12)分别计算k时刻仿射变换后
Figure BDA0002443506160000107
沿
Figure BDA0002443506160000108
的中心坐标τ(k)，
Figure BDA0002443506160000109
沿
Figure BDA00024435061600001010
半轴的平方σ(k)，
Figure BDA00024435061600001011
正交于
Figure BDA00024435061600001012
的半轴的平方δ(k)。
Figure BDA00024435061600001013
Figure BDA00024435061600001014
Figure BDA00024435061600001015
其中ρ为大于零的数，可人为进行设定，一般设置为较小的数值，如ρ＝10-6，
Figure BDA00024435061600001016
为中间变量，
Figure BDA00024435061600001017
为仿射变换后的E，
Figure BDA00024435061600001018
为仿射变换后的φ。
若|μ(k)|≤ρ，则按式(13)至式(15)分别计算k时刻仿射变换后
Figure BDA00024435061600001019
沿
Figure BDA00024435061600001020
的中心坐标τ(k)，
Figure BDA00024435061600001021
沿
Figure BDA00024435061600001022
半轴的平方σ(k)，
Figure BDA00024435061600001023
正交
Figure BDA00024435061600001024
的半轴的平方δ(k)。
τ(k)＝0 (13)
σ(k)＝nα2 (14)
Figure BDA00024435061600001025
其中
Figure BDA0002443506160000111
为中间变量。
S3、按式(16)更新k时刻椭球集合E(k)的中心θc(k)，按式(17)更新k时刻椭球集合E(k)的轴信息矩阵P(k)。
Figure BDA0002443506160000112
Figure BDA0002443506160000113
步骤104，根据弹簧阻尼系统参数向量对应的椭球集合，确定椭球集合对应的支持正多胞体集合和支持正多胞体的交集。
根据k时刻弹簧阻尼系统参数向量对应的椭球集合E(k)，按式(18)确定椭球集合对应的支持多胞体集合O(k)：
Figure BDA0002443506160000114
具体地，
Figure BDA0002443506160000115
Figure BDA0002443506160000116
其中，
Figure BDA0002443506160000117
<·>为内积函数。
根据k时刻椭球集合对应的支持正多胞体集合O(k)，确定支持正多胞体的交集X(k)：
X(k)＝O(1)∩…∩O(k)＝X(k-1)∩O(k) (21)
具体地，按式(22)确定支持正多胞体的交集X(k)。
Figure BDA0002443506160000118
具体地，
Figure BDA0002443506160000119
步骤105，根据弹簧阻尼系统的故障指示信号，确定系统的故障状态和故障时间。
若k时刻故障指示信号f(k)＝1，则表明弹簧阻尼系统在k时刻发生故障；
若k时刻故障指示信号f(k)＝0，则表明弹簧阻尼系统在k时刻未发生故障。
步骤106，若系统发生故障，确定弹簧阻尼系统的参数向量θ的具体故障分量。
步骤1061，对k-1时刻支持向量的交集X(k-1)按式(23)进行n-1维扩展，获得测试集合
Figure BDA0002443506160000121
即对于第i个测试集合
Figure BDA0002443506160000122
除掉第i维，都进行扩展。
Figure BDA0002443506160000123
其中，
Figure BDA0002443506160000124
为第u维参数θu的最大变化范围，u∈{1,…,n}。
步骤1062，根据k-1时刻的测试集合
Figure BDA0002443506160000125
计算k-1时刻测试集合对应的
Figure BDA0002443506160000126
-John椭球
Figure BDA0002443506160000127
具体地：
根据k-1时刻的测试集合
Figure BDA0002443506160000128
的顶点
Figure BDA0002443506160000129
按如下约束求解
Figure BDA00024435061600001210
John椭球
Figure BDA00024435061600001211
的中心
Figure BDA00024435061600001212
Figure BDA00024435061600001213
Qi∈Rn×n为对称正定矩阵，
Figure BDA00024435061600001214
为
Figure BDA00024435061600001215
的第p个顶点。
按式(25)求解
Figure BDA00024435061600001216
-John椭球
Figure BDA00024435061600001217
的轴信息矩阵
Figure BDA00024435061600001218
Figure BDA00024435061600001219
步骤1063，利用弹簧阻尼系统的系统方程和信息向量，按步骤103中更新参数向量对应的椭球集合的方法更新
Figure BDA00024435061600001220
-John椭球。
利用k时刻物块的位移y(k)、信息向量φ(k)，k-1时刻的
Figure BDA00024435061600001221
-John椭球
Figure BDA00024435061600001222
的中心
Figure BDA00024435061600001223
轴信息矩阵Pi t(k-1)和系统不确定噪声的边界γ，按步骤1031和步骤1033的递推更新方法来更新k时刻
Figure BDA00024435061600001224
-John椭球
Figure BDA00024435061600001225
的中心
Figure BDA00024435061600001226
和轴信息矩阵Pi t(k)，i∈{1,…,n}；
利用k+1时刻物块的位移y(k+1)、信息向量φ(k+1)，k时刻的
Figure BDA00024435061600001227
-John椭球
Figure BDA00024435061600001228
的中心
Figure BDA00024435061600001229
轴信息矩阵Pi t(k)和系统不确定噪声的边界γ，按步骤1031和步骤1033的递推更新方法来更新k+1时刻
Figure BDA00024435061600001232
-John椭球
Figure BDA00024435061600001230
的中心
Figure BDA00024435061600001231
和轴信息矩阵Pi t(k+1)，i∈{1,…,n}；
……
利用k+L时刻物块的位移y(k+L)、信息向量φ(k+L)，k+L-1时刻的
Figure BDA00024435061600001328
-John椭球
Figure BDA0002443506160000131
的中心
Figure BDA0002443506160000132
轴信息矩阵Pi t(k+L-1)和系统不确定噪声的边界γ，按步骤1031和步骤1033的递推更新方法来更新k+L时刻
Figure BDA00024435061600001326
-John椭球
Figure BDA0002443506160000133
的中心
Figure BDA0002443506160000134
和轴信息矩阵Pi t(k+L)，i∈{1,…,n}。
L的值是预先设置的，L≤N-k，L的值根据实际需要确定，比如L＝10。
步骤1064，计算k+L时刻
Figure BDA00024435061600001325
-John椭球
Figure BDA0002443506160000135
空集指示信号
Figure BDA0002443506160000136
利用k+L时刻物块的位移y(k+L)、信息向量φ(k+L)，k+L-1时刻的
Figure BDA00024435061600001329
-John椭球
Figure BDA0002443506160000137
的中心
Figure BDA0002443506160000138
轴信息矩阵Pi t(k+L-1)和系统不确定噪声的边界γ，按步骤1031和步骤1032中确定弹簧阻尼系统的故障指示信号的方法来计算k+L时刻
Figure BDA00024435061600001333
-John椭球
Figure BDA0002443506160000139
的空集指示信号
Figure BDA00024435061600001310
步骤1065，根据k+L时刻
Figure BDA00024435061600001334
-John椭球
Figure BDA00024435061600001311
空集指示信号
Figure BDA00024435061600001312
确定弹簧阻尼系统的参数向量θ的具体故障分量。
具体地：
若对于所有的i，i∈{1,…,n}，当k+L时刻
Figure BDA00024435061600001332
-John椭球
Figure BDA00024435061600001313
空集指示信号
Figure BDA00024435061600001314
都为1时，则弹簧阻尼系统的参数向量θ中的θi，i∈{1,…,n}，都发生故障；
若对于所有的i，i∈{1,…,n}，当i≠j，j∈{1,…,n}时，k+L时刻
Figure BDA00024435061600001330
-John椭球
Figure BDA00024435061600001315
空集指示信号
Figure BDA00024435061600001316
都为1，同时当i＝j时，k+L时刻
Figure BDA00024435061600001331
-John椭球
Figure BDA00024435061600001317
空集指示信号
Figure BDA00024435061600001318
为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 BDA00024435061600001335
-John椭球
Figure BDA00024435061600001319
空集指示信号
Figure BDA00024435061600001320
都为1，同时当i≠j且i≠q时，k+L时刻
Figure BDA00024435061600001327
-John椭球
Figure BDA00024435061600001321
空集指示信号
Figure BDA00024435061600001322
都为0，则弹簧阻尼系统的参数向量θ中的θi，i∈{1,…,n}，i＝j且i＝q都发生故障；
若对于所有的i，i∈{1,…,n}，当i＝j，j∈{1,…,n}时，k+L时刻
Figure BDA00024435061600001336
-John椭球
Figure BDA00024435061600001323
空集指示信号
Figure BDA00024435061600001324
为1，同时当i≠j时，k+L时刻
Figure BDA00024435061600001337
-John椭球Ei t(k+L)空集指示信号
Figure BDA00024435061600001411
都为0，则弹簧阻尼系统的参数向量θ中的θi，i∈{1,…,n}，i＝j发生故障。
步骤107，根据弹簧阻尼系统的参数向量的具体故障分量确定扩展方向，按扩展方向重置支持正多胞体的交集。
根据弹簧阻尼系统的参数向量θ的具体故障分量确定扩展方向，按扩展方向重置k-1时刻支持正多胞体的交集。
若弹簧阻尼系统的参数向量θ中的θi，i∈{1,…,n}发生故障，则按式(26)和式(27)分别更新k-1时刻重置的支持正多胞体的交集Xr(k-1)的第i维参数的上界和下界：
Figure BDA0002443506160000141
Figure BDA0002443506160000142
其中，δi max为第i维参数θi的最大变化范围，i∈{1,…,n}，r为上标，带有此上标的参数和集合表示重置后的参数和集合。
若弹簧阻尼系统的参数向量θ中的θi，i∈{1,…,n}未发生故障，则按式(28)和式(29)分别更新k-1时刻重置的支持正多胞体的交集Xr(k-1)的第i维参数的上界和下界：
Figure BDA0002443506160000143
Figure BDA0002443506160000144
根据k-1时刻重置的支持正多胞体的交集Xr(k-1)每维参数的上界和下界，按式(30)获得k-1时刻按扩展方向重置的支持正多胞体的交集Xr(k-1)：
Figure BDA0002443506160000145
步骤108，根据按扩展方向重置的支持正多胞体的交集，识别弹簧阻尼系统的故障参数。
步骤1081，根据k-1时刻按扩展方向重置的支持正多胞体的交集Xr(k-1)，按步骤1062中计算k-1时刻测试集合
Figure BDA0002443506160000146
对应的
Figure BDA0002443506160000148
-John椭球
Figure BDA0002443506160000147
的方法计算k-1时刻按扩展方向重置的支持正多胞体的交集Xr(k-1)对应的
Figure BDA0002443506160000149
-John椭球Er(k-1)。
步骤1082，利用k时刻物块的位移y(k)、信息向量φ(k)，k-1时刻的
Figure BDA00024435061600001410
-John椭球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 BDA0002443506160000151
Figure BDA0002443506160000152
Figure BDA0002443506160000153
为k时刻支持正多胞体的交集X(k)的第i维对应的参数的上界，
Figure BDA0002443506160000154
为k时刻支持正多胞体的交集X(k)的第i维对应的参数的下界，i∈{1,…,n}。
利用k时刻支持正多胞体的交集X(k)的每一维对应的参数的中心值
Figure BDA0002443506160000155
按式(32)确定故障参数向量
Figure BDA0002443506160000156
从而实现k时刻已发生故障的弹簧阻尼系统的故障参数的识别。
Figure BDA0002443506160000157
需要说明的是：本发明实施例提供的面向工业生产过程的弹簧阻尼系统滤波故障诊断方法，诊断弹簧阻尼系统是否发生故障时在弹簧阻尼系统处于工作状态下进行的。
为验证本申请所提出的面向工业生产过程的弹簧阻尼系统滤波故障诊断方法的收敛速度和故障识别时间，进行如下仿真实验：
仿真实验中设置系统噪声为|e(k)|≤0.01，设定的弹簧阻尼系统的相关状态和参数如表1所示，即当时间k到达1001，2001和3001时，我们分别加入故障状态1，故障状态2和故障状态3。
表1弹簧阻尼系统在不同状态下的参数分量的数值
Figure BDA0002443506160000158
Figure BDA0002443506160000161
基于相同的仿真条件，将本申请提出的一种面向工业生产过程的弹簧阻尼系统滤波故障诊断方法与基于全域放大集员滤波的弹簧阻尼系统故障诊断方法进行对比，得到的故障诊断对比结果如图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 BDA0002443506160000171
需要说明的是，表2中
Figure BDA0002443506160000172
表示向工业生产过程的弹簧阻尼系统滤波故障诊断方法所得出的参数识别结果，
Figure BDA0002443506160000173
表示基于全域放大集员滤波的弹簧阻尼系统故障诊断方法所得出的参数识别结果。对比表1可知，这两种算法在识别精度上都比较高。但是结合图3-图6可知，本申请所提出的工业生产过程的弹簧阻尼系统滤波故障诊断方法能更快地识别出故障。
综上所述，本申请针对工业生产过程中弹簧阻尼系统故障诊断问题，提出了一类面向工业生产过程的弹簧阻尼系统滤波故障诊断方法，通过建立弹簧阻尼系统的离散系统模型，确定弹簧阻尼系统的系统方程；获取弹簧阻尼系统在工作状态下的外加控制力和物块的位移，确定系统的信息向量；根据弹簧阻尼系统的系统方程和信息向量，确定参数向量对应的椭球集合和系统故障指示信号的数值；根据弹簧阻尼系统参数向量对应的椭球集合，确定椭球集合对应的支持正多胞体集合和支持正多胞体的交集；根据弹簧阻尼系统的故障指示信号，确定系统的故障状态和故障时间；若系统发生故障，确定弹簧阻尼系统的参数向量的具体故障分量；根据弹簧阻尼系统的参数向量的具体故障分量，按扩展方向重置支持正多胞体的交集；根据按扩展方向重置的支持正多胞体的交集，识别弹簧阻尼系统的故障参数；解决了系统噪声不确定下的弹簧阻尼系统的快速准确实时的故障诊断，同时相比于基于全域放大集员滤波的故障诊断方法，面向工业生产过程的弹簧阻尼系统滤波故障诊断方法在弹簧阻尼系统的故障诊断中有着收敛速度快，故障识别时间短的优点。
本发明实施例中的部分步骤，可以利用软件实现，相应的软件程序可以存储在可读取的存储介质中，如光盘或硬盘等。
以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种弹簧阻尼系统滤波故障诊断方法，其特征在于，所述方法包括：
根据弹簧阻尼系统的系统方程和信息向量确定弹簧阻尼系统的参数向量对应的椭球集合，进而确定椭球集合对应的支持正多胞体集合和支持正多胞体的交集；
在系统发生故障时，根据测试集合及其
Figure FDA0002443506150000011
-John椭球确定弹簧阻尼系统的参数向量的具体故障分量，继而根据弹簧阻尼系统的参数向量的具体故障分量确定扩展方向，按扩展方向重置支持正多胞体的交集；
根据按扩展方向重置的支持正多胞体的交集，识别弹簧阻尼系统的参数向量中的故障参数。
2.根据权利要求1所述的方法，其特征在于，所述方法包括：
步骤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中按扩展方向重置的支持正多胞体的交集，识别弹簧阻尼系统的故障参数。
3.根据权利要求2所述的方法，其特征在于，所述步骤104，根据步骤103确定的弹簧阻尼系统参数向量θ对应的椭球集合，确定椭球集合对应的支持正多胞体集合O(k)和支持正多胞体的交集，包括：
Figure FDA0002443506150000021
Figure FDA0002443506150000022
Figure FDA0002443506150000023
其中，
Figure FDA0002443506150000024
表示k时刻支持正多胞体集合O(k)的参数的上界，
Figure FDA0002443506150000025
表示k时刻支持正多胞体集合O(k)的参数的下界，
Figure FDA0002443506150000026
u∈{1,2,…,n}，<·>为内积函数，n为参数向量θ的维数，k时刻表示第k个离散时刻；
根据支持正多胞体集合O(k)，确定k时刻支持正多胞体的交集：
X(k)＝O(1)∩…∩O(k)＝X(k-1)∩O(k)
即
Figure FDA0002443506150000027
Figure FDA0002443506150000028
Figure FDA0002443506150000029
其中，
Figure FDA00024435061500000210
表示k时刻支持正多胞体的交集X(k)的参数的上界，
Figure FDA00024435061500000211
表示k时刻支持正多胞体的交集X(k)的参数的下界。
4.根据权利要求3所述的方法，其特征在于，所述步骤106，若系统发生故障，确定弹簧阻尼系统的参数向量θ的具体故障分量，包括：
对k-1时刻支持正多胞体的交集X(k-1)进行n-1维扩展，获得测试集合
Figure FDA00024435061500000212
i∈{1,…,n}，即对于第i个测试集合
Figure FDA00024435061500000213
除掉第i维，都进行扩展；
根据所述k-1时刻的测试集合
Figure FDA00024435061500000214
计算k-1时刻测试集合对应的
Figure FDA00024435061500000215
-John椭球
Figure FDA00024435061500000216
利用弹簧阻尼系统的系统方程和信息向量，更新k-1时刻的
Figure FDA00024435061500000222
-John椭球
Figure FDA00024435061500000217
得到k时刻的
Figure FDA00024435061500000218
-John椭球
Figure FDA00024435061500000219
计算所述k+L时刻
Figure FDA00024435061500000223
-John椭球
Figure FDA00024435061500000220
空集指示信号
Figure FDA00024435061500000221
根据所述k+L时刻
Figure FDA0002443506150000031
-John椭球
Figure FDA0002443506150000032
空集指示信号
Figure FDA0002443506150000033
确定弹簧阻尼系统的参数向量θ的具体故障分量。
5.根据权利要求4所述的方法，其特征在于，所述步骤107，根据步骤106确定的弹簧阻尼系统的参数向量的具体故障分量确定扩展方向，按扩展方向重置支持正多胞体的交集，包括：
若弹簧阻尼系统的参数向量θ中的θi发生故障，按下式更新k-1时刻重置的支持正多胞体的交集Xr(k-1)的第i维参数的上界：
Figure FDA0002443506150000034
若弹簧阻尼系统的参数向量θ中的θi发生故障，按下式更新k-1时刻重置的支持正多胞体的交集Xr(k-1)的第i维参数的下界：
Figure FDA0002443506150000035
若弹簧阻尼系统的参数向量θ中的θi未发生故障，按下式更新k-1时刻重置的支持正多胞体的交集Xr(k-1)的第i维参数的上界：
Figure FDA0002443506150000036
若弹簧阻尼系统的参数向量θ中的θi未发生故障，按下式更新k-1时刻重置的支持正多胞体的交集Xr(k-1)的第i维参数的下界：
Figure FDA0002443506150000037
根据所述k-1时刻支持正多胞体的交集Xr(k-1)每维参数的上界和下界，按下式获得k-1时刻按扩展方向重置的支持正多胞体的交集Xr(k-1)：
Figure FDA0002443506150000038
其中，
Figure FDA0002443506150000039
为第i维参数θi的最大变化范围，i∈{1,…,n}，r为上标，带有此上标的参数和集合表示重置后的参数和集合。
6.根据权利要求5所述的方法，其特征在于，所述根据步骤107中按扩展方向重置的支持正多胞体的交集，识别弹簧阻尼系统的故障参数，包括：
根据所述k-1时刻按扩展方向重置的支持正多胞体的交集Xr(k-1)，计算k-1时刻按扩展方向重置的支持正多胞体的交集Xr(k-1)对应的
Figure FDA0002443506150000041
-John椭球Er(k-1)；
利用k时刻物块的位移y(k)、信息向量φ(k)，k-1时刻的
Figure FDA00024435061500000415
-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 FDA0002443506150000042
Figure FDA0002443506150000043
其中
Figure FDA0002443506150000044
i∈{1,…,n}，
Figure FDA0002443506150000045
为k时刻支持正多胞体的交集X(k)的第i维对应的参数的上界，
Figure FDA0002443506150000046
为k时刻支持正多胞体的交集X(k)的第i维对应的参数的下界。
7.根据权利要求6所述的方法，其特征在于，所述步骤103，根据步骤101确定的弹簧阻尼系统的系统方程和步骤102确定的信息向量φ(k)，确定参数向量θ对应的椭球集合和系统故障指示信号f(k)的数值，包括：
按下面两式确定k时刻仿射变换后
Figure FDA0002443506150000047
正交的两个平行超平面的
Figure FDA0002443506150000048
坐标，即α(k)和
Figure FDA0002443506150000049
其中
Figure FDA00024435061500000410
为仿射变换后的φ：
Figure FDA00024435061500000411
Figure FDA00024435061500000412
若α(k)≥1或
Figure FDA00024435061500000413
则故障指示信号f(k)＝1，表明弹簧阻尼系统在k时刻发生故障；
θc(k)＝θc(k-1)，
P(k)＝P(k-1)；
若α(k)≤1且
Figure FDA00024435061500000414
则故障指示信号f(k)＝0，表明弹簧阻尼系统在k时刻未发生故障。
8.根据权利要求7所述的方法，其特征在于，若
Figure FDA0002443506150000051
按下面两式更新k时刻椭球集合E(k)的中心θc(k)和轴信息矩阵P(k)：
θc(k)＝θc(k-1)，
P(k)＝P(k-1)；
若两者还同时满足
Figure FDA0002443506150000052
则
Figure FDA0002443506150000053
在上述条件下，若|μ(k)|＞ρ，则
Figure FDA0002443506150000054
Figure FDA0002443506150000055
Figure FDA0002443506150000056
Figure FDA0002443506150000057
Figure FDA0002443506150000058
在上述条件下，若|μ(k)|≤ρ，则
Figure FDA0002443506150000059
τ(k)＝0，
σ(k)＝nα2，
Figure FDA00024435061500000510
之后按下面两式更新k时刻椭球集合E(k)的中心θc(k)和轴信息矩阵P(k)
Figure FDA00024435061500000511
Figure FDA00024435061500000512
其中n为参数向量θ的维数，μ(k)为α(k)和
Figure FDA00024435061500000513
的平均值，ρ为大于零的数，设定ρ＝10-6，τ(k)为k时刻仿射变换后
Figure FDA0002443506150000061
沿
Figure FDA0002443506150000062
的中心坐标，σ(k)为
Figure FDA0002443506150000063
沿
Figure FDA0002443506150000064
半轴的平方，δ(k)为
Figure FDA0002443506150000065
正交于
Figure FDA0002443506150000066
的半轴的平方，ε(k)、b(k)、α为中间变量，
Figure FDA0002443506150000067
为仿射变换后的E，
Figure FDA0002443506150000068
为仿射变换后的φ。
9.根据权利要求8所述的方法，其特征在于，所述步骤102，获取弹簧阻尼系统在工作状态下的外加控制力和对应的物块的位移，以确定系统的信息向量φ(k)，包括：
在预定时间范围内，获取弹簧阻尼系统在工作状态下的外加控制力和对应的物块的位移；
将所得的外加控制力和对应的物块的位移的数据代入下式中：
φ(k)＝[-y(k-1),-y(k-2),u(k-1),u(k-2)]T
确定出弹簧阻尼系统的信息向量φ(k)；k的取值范围为1至N，k为整数。
10.一种弹簧阻尼系统滤波故障诊断系统，其特征在于，所述系统采用权利要求1-9任一所述的方法对弹簧阻尼系统进行故障诊断。
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 true CN111597647A (zh) 2020-08-28
CN111597647B 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)

Cited By (4)

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

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" *

Cited By (8)

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

Also Published As

Publication number Publication date
CN111597647B (zh) 2023-04-25

Similar Documents

Publication Publication Date Title
CN111597647A (zh) 一种面向工业生产过程的弹簧阻尼系统滤波故障诊断方法
CN111639678B (zh) 一种基于集成神经网络的ins/gps组合导航故障检测与诊断方法
CN107656265B (zh) 针对多帧检测前跟踪短航迹的粒子滤波融合方法
CN112001110B (zh) 一种基于振动信号空间时时递归图卷积神经网络的结构损伤识别监测方法
CN113111820B (zh) 基于改进cnn与关系模块的旋转部件故障诊断方法及其装置
Camelio et al. Sensor placement for effective diagnosis of multiple faults in fixturing of compliant parts
CN113218537B (zh) 温度异常检测模型的训练方法、装置、设备和存储介质
CN109623489B (zh) 一种改进的机床健康状态评定方法及数控机床
CN112926152B (zh) 一种数字孪生驱动的薄壁件装夹力精准控制与优化方法
CN111680725A (zh) 基于重构贡献的气体传感器阵列多故障隔离算法
Wang et al. Fault detection based on diffusion maps and k nearest neighbor diffusion distance of feature space
CN115290316A (zh) 一种偏心旋转阀故障诊断方法
Liu et al. A novel strategy for response and force reconstruction under impact excitation
CN100465997C (zh) 基于元胞自动机的图像边缘检测算法
CN116909901B (zh) 基于富样本径向投影可信聚类的软件性能可视评价方法
Li et al. Symbolization‐based differential evolution strategy for identification of structural parameters
CN116630320A (zh) 电池极片的检测方法、装置、电子设备及存储介质
CN109543626B (zh) 一种机械装备旋转部件轴心轨迹识别方法
CN112069998A (zh) 基于卷积神经网络的gbdt故障诊断集成方法
CN116451470A (zh) 基于多元线性回归和动力学的机器人能耗建模方法及装置
CN116881781A (zh) 一种运行模态阻尼识别方法、损伤检测方法、系统及设备
CN116486146A (zh) 一种旋转机械设备的故障检测方法、系统、装置和介质
CN112747773B (zh) 基于Allan方差和随机多项式提高陀螺仪精度的方法
CN115114985A (zh) 一种基于集合理论的传感器系统分布式融合方法
CN104516858A (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