CN113236506B - 一种基于滤波的工业时滞系统故障检测方法 - Google Patents

一种基于滤波的工业时滞系统故障检测方法 Download PDF

Info

Publication number
CN113236506B
CN113236506B CN202110550989.3A CN202110550989A CN113236506B CN 113236506 B CN113236506 B CN 113236506B CN 202110550989 A CN202110550989 A CN 202110550989A CN 113236506 B CN113236506 B CN 113236506B
Authority
CN
China
Prior art keywords
time
lag
fully
industrial
matrix
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
CN202110550989.3A
Other languages
English (en)
Other versions
CN113236506A (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 CN202110550989.3A priority Critical patent/CN113236506B/zh
Publication of CN113236506A publication Critical patent/CN113236506A/zh
Application granted granted Critical
Publication of CN113236506B publication Critical patent/CN113236506B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F03—MACHINES OR ENGINES FOR LIQUIDS; WIND, SPRING, OR WEIGHT MOTORS; PRODUCING MECHANICAL POWER OR A REACTIVE PROPULSIVE THRUST, NOT OTHERWISE PROVIDED FOR
    • F03D—WIND MOTORS
    • F03D17/00—Monitoring or testing of wind motors, e.g. diagnostics

Landscapes

  • Engineering & Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Sustainable Development (AREA)
  • Sustainable Energy (AREA)
  • Chemical & Material Sciences (AREA)
  • Combustion & Propulsion (AREA)
  • Mechanical Engineering (AREA)
  • General Engineering & Computer Science (AREA)
  • Control Of Eletrric Generators (AREA)

Abstract

本发明公开了一种基于滤波的工业时滞系统故障检测方法，属于时滞系统故障检测技术领域。所述方法根据所建立的工业时滞系统的离散线性模型和获取的系统输出，设计最优全对称多胞体卡尔曼滤波器；再根据所设计的最优全对称多胞体卡尔曼滤波器，估计工业时滞系统的状态向量对应的全对称多胞体；之后根据状态向量对应的全对称多胞体，计算工业时滞系统的输出向量对应的全对称多胞体的上下界；最后根据输出向量对应的上下界，确定工业时滞系统的故障情况，针对时滞系统，本申请通过在设计对应的最优全对称多胞体卡尔曼滤波器时考虑了时滞项的作用，从而针对时滞系统能够得到更准确的状态估计和故障检测结果。

Description

一种基于滤波的工业时滞系统故障检测方法
技术领域
本发明涉及一种基于滤波的工业时滞系统故障检测方法，属于时滞系统故障检测技术领域。
背景技术
工业系统一般指工业控制系统，是由各种自动化控制组件以及对实时数据进行采集、监测的过程控制组件共同构成的确保工业基础设施自动化运行、过程控制与监控的业务流程管控系统。各类工业系统工作环境复杂，受各种不可控因素影响，系统易发生故障，因此对系统进行可行有效的故障检测，保证系统安全可靠的运行时非常有必要的。
实际应用中，随着系统日趋复杂化和规模化，系统噪声干扰难以满足特定的分布规律的要求，因此采用有效的基于集员滤波的故障检测方法对噪声未知但有界的系统进行故障检测是有着重要意义的。
但同时，在实际工业生产过程中，受信息传输等因素的影响，各类工业系统不可避免地会出现时滞的现象，给系统的稳定性能带来不利影响，同时也影响着系统状态估计和故障检测的结果。因此，如何减小系统时滞对状态估计和故障检测带来的影响，保证工业时滞系统的故障检测精确度是当前故障检测领域亟待解决的重要问题。
现有技术中，对于噪声未知但有界的工业系统，通常能够采用合适的基于集员滤波的故障检测方法对系统的故障进行检测，但是针对更复杂的时滞系统的故障检测，没有考虑时滞作用的算法在进行状态估计和故障检测过程中就忽略了时滞项，在迭代运算中不考虑此时滞项起到的作用，因此得到的结果肯定与原有的时滞系统肯定有出入，所以系统将会出现状态估计和故障检测结果不准确的情况。
发明内容
为了实现对工业时滞系统故障的准确检测，本发明提供了一种基于滤波的工业时滞系统故障检测方法，所述方法包括：
步骤101，根据工业时滞系统的工作原理，建立工业时滞系统的离散线性模型；
步骤102，获取工业时滞系统在工作状态下的输出向量；
步骤103，根据所建立的工业时滞系统的离散线性模型和工业时滞系统在工作状态下的输出向量，设计针对工业时滞系统的最优全对称多胞体卡尔曼滤波器；
步骤104，根据所设计的针对工业时滞系统的最优全对称多胞体卡尔曼滤波器，估计工业时滞系统的状态向量对应的全对称多胞体；
步骤105，根据状态向量对应的全对称多胞体，计算工业时滞系统的输出向量对应的全对称多胞体的上下界；
步骤106，根据输出向量对应的上下界，确定工业时滞系统的故障情况。
可选的，所述步骤101建立的建立工业时滞系统的离散线性模型为：
Figure BDA0003072478400000021
其中x(k)，u(k)和y(k)分别表示工业时滞系统在k时刻的状态向量，输入向量和输出向量，k为离散时间；
w(k)和v(k)分别为过程噪声和测量噪声，并且均有界，分别为w(k)∈W和v(k)∈V，其中W为过程噪声w(k)对应的全对称多胞体，V为测量噪声v(k)对应的全对称多胞体；
A，B，C，D，F分别为对应的参数矩阵，Ah为时滞项参数矩阵，h为时滞时间常数。
可选的，所述步骤103包括：
当0≤k≤h时，工业时滞系统的系统时滞对状态估计结果没有影响，按式(2)～(8)设计最优全对称多胞体卡尔曼滤波器：
Figure BDA0003072478400000022
Figure BDA0003072478400000023
Figure BDA0003072478400000024
Figure BDA0003072478400000025
K(k)＝R(k)S-1(k) (6)
L(k)＝AK(k) (7)
Figure BDA0003072478400000026
其中，
Figure BDA00030724784000000212
为
Figure BDA0003072478400000029
的降维矩阵，
Figure BDA00030724784000000210
为k时刻状态向量对应的全对称多胞体的生成矩阵，L(k)为k时刻的最优增益矩阵，Gv为测量噪声v(k)对应的全对称多胞体V的生成矩阵，Qv、
Figure BDA00030724784000000211
S(k)、R(k)、K(k)为中间计算变量，“^”表示估计值符号；
当k＞h时，工业时滞系统的系统时滞对状态估计结果存在影响，按式(9)～(20)设计最优全对称多胞体卡尔曼滤波器：
Figure BDA0003072478400000031
Figure BDA0003072478400000032
Figure BDA0003072478400000033
Figure BDA0003072478400000034
Figure BDA0003072478400000035
Figure BDA0003072478400000036
Figure BDA0003072478400000037
Figure BDA0003072478400000038
K1(k)＝R1(k)S-1(k) (17)
K2(k)＝R2(k)S-1(k) (18)
L(k)＝AK1(k)+AhK2(k) (19)
Figure BDA0003072478400000039
其中“∑”表示累加符号，“∏”表示累乘符号，Gw为w(k)对应的全对称多胞体W的生成矩阵，
Figure BDA00030724784000000310
为
Figure BDA00030724784000000311
的降维矩阵，
Figure BDA00030724784000000312
为k-h时刻状态向量对应的全对称多胞体的生成矩阵，
Figure BDA00030724784000000313
为
Figure BDA00030724784000000314
的降维矩阵，
Figure BDA00030724784000000315
为k-h-i时刻状态向量对应的全对称多胞体的生成矩阵，L(k-i)为k-i时刻的最优增益矩阵，
Figure BDA00030724784000000316
R1(k)、R2(k)、K1(k)、K2(k)为中间计算变量。
可选的，所述步骤104包括：
当0≤k≤h时，状态向量对应的全对称多胞体
Figure BDA0003072478400000041
具体可通过式(21)和(22)计算得到：
Figure BDA0003072478400000042
Figure BDA0003072478400000043
其中
Figure BDA0003072478400000044
为k+1时刻状态向量对应的全对称多胞体，
Figure BDA0003072478400000045
和
Figure BDA0003072478400000046
为别为
Figure BDA0003072478400000047
的中心和生成矩阵，
Figure BDA0003072478400000048
和
Figure BDA0003072478400000049
分别为k-h时刻状态向量对应的全对称多胞体的中心和生成矩阵，且
Figure BDA00030724784000000410
当k＞h时，状态向量对应的全对称多胞体
Figure BDA00030724784000000411
具体可通过式(23)～(30)计算得到：
Figure BDA00030724784000000412
Figure BDA00030724784000000413
Figure BDA00030724784000000414
Figure BDA00030724784000000415
Figure BDA00030724784000000416
Figure BDA00030724784000000417
Figure BDA00030724784000000418
Figure BDA00030724784000000419
其中H1(k)、H1i(k)、H2(k)、H2i(k)、H3(k)、H3i(k)为中间计算变量，为状态向量x(k)的维度，r为全对称多胞体
Figure BDA00030724784000000420
的生成矩阵的维度，
Figure BDA00030724784000000421
为H1(k)的维度，
Figure BDA00030724784000000422
为H2(k)的维度，
Figure BDA0003072478400000051
为H3(k)的维度。
可选的，所述步骤105包括：
根据状态向量对应的全对称多胞形
Figure BDA0003072478400000052
基于测量噪声v(k)对应的全对称多胞体V的生成矩阵Gv，按式(31)计算k时刻输出矩阵对应的全对称多胞体
Figure BDA0003072478400000053
Figure BDA0003072478400000054
其中
Figure BDA0003072478400000055
为k时刻锂电池系统输出矩阵对应的全对称多胞体
Figure BDA0003072478400000056
的中心，
Figure BDA0003072478400000057
为k时刻锂电池系统输出矩阵对应的全对称多胞体
Figure BDA0003072478400000058
的生成矩阵，其中
Figure BDA0003072478400000059
表示闵可夫斯基和，“⊙”表示线性映射；
根据输出矩阵对应的全对称多胞体
Figure BDA00030724784000000510
按式(32)和(33)计算k时刻工业时滞系统输出向量对应的上下界：
Figure BDA00030724784000000511
Figure BDA00030724784000000512
其中
Figure BDA00030724784000000513
表述输出向量对应的上界，
Figure BDA00030724784000000514
表示输出向量对应的下界。
可选的，所述步骤106包括：
若步骤105得出的输出向量的上下界满足
Figure BDA00030724784000000515
和
Figure BDA00030724784000000516
则表明工业时滞系统未发生故障，此时故障检测信号f(k)＝0；
否则，表明该工业时滞系统在k时刻发生故障，此时故障检测信号f(k)＝1；
其中
Figure BDA00030724784000000517
和
Figure BDA00030724784000000518
分别为输出矩阵
Figure BDA00030724784000000519
和
Figure BDA00030724784000000520
中的元素。
可选的，所述工业时滞系统包括时滞风电机组桨距子系统。
可选的，所述方法应用于时滞风电机组桨距子系统时，所述步骤101之前还包括：
根据风电机组桨距子系统的工作原理，建立时滞风电机组桨距子系统的系统模型：
Figure BDA0003072478400000061
其中，β为桨距角，βa为桨速度，βr表示桨距参考值，ωn和ζ为时滞风电机组桨距子系统的自然频率和阻尼系数，“·”表示求导符号。
可选的，所述时滞风电机组桨距子系统表示为连续时间状态空间方程为：
Figure BDA0003072478400000062
其中x＝[β,βα]T为时滞风电机组桨距子系统的状态向量，u＝βr为时滞风电机组桨距子系统的输入向量，w和v分别表示时滞风电机组桨距子系统的系统的过程噪声和测量噪声。
可选的，针对时滞风电机组桨距子系统，
Figure BDA0003072478400000063
Figure BDA0003072478400000064
为时滞项参数矩阵，h＝6为时滞时间常数。
本发明有益效果是：
本申请根据所建立的工业时滞系统的离散线性模型和获取的系统输出，设计最优全对称多胞体卡尔曼滤波器；再根据所设计的最优全对称多胞体卡尔曼滤波器，估计工业时滞系统的状态向量对应的全对称多胞体；之后根据状态向量对应的全对称多胞体，计算工业时滞系统的输出向量对应的全对称多胞体的上下界；最后根据输出向量对应的上下界，确定工业时滞系统的故障情况，针对时滞系统，本申请通过在设计对应的最优全对称多胞体卡尔曼滤波器时考虑了时滞项的作用，从而针对时滞系统能够得到更准确的状态估计和故障检测结果。
附图说明
为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1是本发明一个实施例公开的一种基于滤波的工业时滞系统故障检测方法流程图。
图2是采用本申请提供的基于滤波的工业时滞系统故障检测方法对时滞风电机组桨距子系统进行故障估计时输出的上下界和实际输出曲线图。
图3是采用本申请提供的基于滤波的工业时滞系统故障检测方法对时滞风电机组桨距子系统进行故障检测的仿真曲线图。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地详细描述。
实施例一：
本实施例提供一种基于滤波的工业时滞系统故障检测方法，参见图1，所述方法包括：
步骤101，根据工业时滞系统的工作原理，建立工业时滞系统的离散线性模型。
首先根据工业时滞系统的工作原理，建立工业时滞系统的系统模型。
设置合适的采样时间Ts对系统进行离散化，同时考虑系统时滞作用，建立工业时滞系统的离散线性模型：
Figure BDA0003072478400000071
其中x(k)，u(k)和y(k)分别表示工业时滞系统在k时刻的状态向量、输入向量和输出向量，k为离散时间；
w(k)和v(k)分别为过程噪声和测量噪声，并且均有界，具体为w(k)∈W和v(k)∈V，其中W为过程噪声w(k)对应的全对称多胞体，V为测量噪声v(k)对应的全对称多胞体；
A，B，C，D，F分别为适当维度的参数矩阵，Ah为时滞项参数矩阵，h为时滞时间常数。
系统初始状态满足x(k-h)∈X(k-h)＝＜p(k-h),G(k-h)＞， 0≤k≤h，其中X(k-h)为系统在k-h时刻的状态向量x(k-h)对应的全对称多胞体，p(k-h)和G(k-h)分别为X(k-h)的中心和生成矩阵。
步骤102，获取工业时滞系统在工作状态下的输出向量。
在预定时间范围内，获取工业时滞系统在工作状态下的测量数据。
预定时间范围为1至N，N为整数，N的值是预先设置的。
实际应用中，可选用合适的传感器或测量仪器获取系统测量数据。
将所测得的工作状态下的系统测量数据代入式(1)所示的离散线性模型中输出向量y(k)中，确定工作状态下工业时滞系统的输出向量y(k)，k的取值范围为1至N，k为整数。
步骤103，根据步骤101所建立的工业时滞系统的离散线性模型和步骤102所获取的工作状态下工业时滞系统的输出向量y(k)，设计最优全对称多胞体卡尔曼滤波器。
当0≤k≤h时，系统时滞对状态估计结果没有影响，按式(2)～(8)设计最优全对称多胞体卡尔曼滤波器：
Figure BDA0003072478400000081
Figure BDA0003072478400000082
Figure BDA0003072478400000083
Figure BDA0003072478400000084
K(k)＝R(k)S-1(k) (6)
L(k)＝AK(k) (7)
Figure BDA0003072478400000085
其中，
Figure BDA0003072478400000086
为
Figure BDA0003072478400000087
的降维矩阵，
Figure BDA0003072478400000088
为k时刻状态向量对应的全对称多胞体的生成矩阵，L(k)为k时刻的最优增益矩阵，Gv为测量噪声v(k)对应的全对称多胞体V的生成矩阵，Qv、
Figure BDA0003072478400000089
S(k)、R(k)、K(k)为中间计算变量，“^”表示估计值符号。
当k＞h时，系统时滞对状态估计结果存在作用，按式(9)～(20)设计最优全对称多胞体卡尔曼滤波器：
Figure BDA00030724784000000810
Figure BDA00030724784000000811
Figure BDA00030724784000000812
Figure BDA00030724784000000813
Figure BDA00030724784000000814
Figure BDA00030724784000000815
Figure BDA0003072478400000091
Figure BDA0003072478400000092
K1(k)＝R1(k)S-1(k) (17)
K2(k)＝R2(k)S-1(k) (18)
L(k)＝AK1(k)+AhK2(k) (19)
Figure BDA0003072478400000093
其中“∑”表示累加符号，“Π”表示累乘符号，Gw为w(k)对应的全对称多胞体W的生成矩阵，
Figure BDA0003072478400000094
为
Figure BDA0003072478400000095
的降维矩阵，
Figure BDA0003072478400000096
为k-h时刻状态向量对应的全对称多胞体的生成矩阵，
Figure BDA0003072478400000097
为
Figure BDA0003072478400000098
的降维矩阵，
Figure BDA0003072478400000099
为k-h-i时刻状态向量对应的全对称多胞体的生成矩阵，L(k-i)为k-i时刻的最优增益矩阵，
Figure BDA00030724784000000910
R1(k)、R2(k)、K1(k)、K2(k)为中间计算变量。
步骤104，根据所设计的最优全对称多胞体卡尔曼滤波器，估计工业时滞系统的状态向量对应的全对称多胞体。
当0≤k≤h时，状态向量对应的全对称多胞体
Figure BDA00030724784000000911
具体可通过式(21)和(22)计算得到：
Figure BDA00030724784000000912
Figure BDA00030724784000000913
其中
Figure BDA00030724784000000914
为k+1时刻状态向量对应的全对称多胞体，
Figure BDA00030724784000000915
和
Figure BDA00030724784000000916
为别为
Figure BDA00030724784000000917
的中心和生成矩阵，
Figure BDA00030724784000000918
和
Figure BDA00030724784000000919
分别为k-h时刻状态向量对应的全对称多胞体的中心和生成矩阵，且
Figure BDA00030724784000000920
当k＞h时，状态向量对应的全对称多胞体
Figure BDA00030724784000000921
具体可通过式(23)～(30)计算得到：
Figure BDA0003072478400000101
Figure BDA0003072478400000102
Figure BDA0003072478400000103
Figure BDA0003072478400000104
Figure BDA0003072478400000105
Figure BDA0003072478400000106
Figure BDA0003072478400000107
Figure BDA0003072478400000108
其中H1(k)、H1i(k)、H2(k)、H2i(k)、H3(k)、H3i(k)为中间计算变量，nx、r、
Figure BDA0003072478400000109
Figure BDA00030724784000001010
和
Figure BDA00030724784000001011
分别为状态向量、全对称多胞体
Figure BDA00030724784000001012
的生成矩阵、H1(k)、H2(k)和H3(k)的维度。
步骤105，根据状态向量对应的全对称多胞体，计算工业时滞系统的输出向量对应的全对称多胞体的上下界。
根据状态向量对应的全对称多胞形
Figure BDA00030724784000001013
基于测量噪声v(k)对应的全对称多胞体V的生成矩阵Gv，按式(31)计算k时刻输出矩阵对应的全对称多胞体
Figure BDA00030724784000001014
Figure BDA00030724784000001015
其中
Figure BDA0003072478400000111
为k时刻锂电池系统输出矩阵对应的全对称多胞体
Figure BDA0003072478400000112
的中心，
Figure BDA0003072478400000113
为k时刻锂电池系统输出矩阵对应的全对称多胞体
Figure BDA0003072478400000114
的生成矩阵，其中
Figure BDA0003072478400000115
表示闵可夫斯基和，“⊙”表示线性映射。
根据输出矩阵对应的全对称多胞体
Figure BDA0003072478400000116
按式(32)和(33)计算k时刻工业时滞系统输出向量对应的上下界：
Figure BDA0003072478400000117
Figure BDA0003072478400000118
其中
Figure BDA0003072478400000119
表述输出向量对应的上界，
Figure BDA00030724784000001110
表示输出向量对应的下界。
步骤106，根据输出向量对应的上下界，确定工业时滞系统的故障情况。
若输出向量的上下界满足
Figure BDA00030724784000001111
和
Figure BDA00030724784000001112
表明工业时滞系统未发生故障，此时故障检测信号f(k)＝0；否则，表明该工业时滞系统在k时刻发生故障，此时故障检测信号f(k)＝1。
其中
Figure BDA00030724784000001113
和
Figure BDA00030724784000001114
分别为输出矩阵
Figure BDA00030724784000001115
和
Figure BDA00030724784000001116
中的元素。
实施例二：
本实施例提供一种基于滤波的工业时滞系统故障检测方法，以应用于时滞风电机组桨距子系统为例进行仿真实验如下，以验证本申请所提出的基于滤波的工业时滞系统故障检测方法的正确性和可行性：
根据风电机组桨距子系统的工作原理，建立桨距子系统的系统模型：
Figure BDA00030724784000001117
其中，β为桨距角，βa为桨速度，βr表示桨距参考值，ωn和ζ为已知的系统参数，分别表示系统的自然频率和阻尼系数，“·”表示求导符号。
根据桨距子系统的系统模型，将桨距子系统表示为连续时间状态空间方程为：
Figure BDA0003072478400000121
其中x＝[β,βα]T为状态向量，u＝βr为输入向量，w和v分别表示系统的过程噪声和测量噪声。
设置采样时间Ts＝0.01s对系统进行离散化，同时考虑系统时滞作用，可建立时滞风电机组桨距子系统线性离散模型为：
Figure BDA0003072478400000122
其中x(k)，u(k)和y(k)分别表示桨距子系统在k时刻的状态向量，输入向量和输出向量，k为离散时间；
w(k)和v(k)分别为过程噪声和测量噪声，并且均有界，具体为
Figure BDA0003072478400000123
其中W为w(k)对应的全对称多胞体，V为v(k)对应的全对称多胞体；
A，B，C，D，F分别为适当维度的参数矩阵，当系统处于正常状态时，取ωn＝11.11rad/s，ζ＝0.6，可得离散化后的各参数矩阵具体为：
Figure BDA0003072478400000124
Figure BDA0003072478400000125
为时滞项参数矩阵，h＝6为时滞时间常数。
系统初始状态满足x(k-h)∈X(k-h)＝＜p(k-h)，
Figure BDA0003072478400000126
Figure BDA0003072478400000127
其中X(k-h)为x(k-h)对应的全对称多胞体，p(k-h)和G(k-h)分别为X(k-h)的中心和生成矩阵。
当时滞风电机组桨距子系统处于工作状态下，给定系统输入u(k)＝βr(k)＝1.5sin(6k)+7；
实际应用中，可利用桨距角测量仪和角速度测量器测量桨距子系统中的桨距角和角速度的大小。
将所得的工作状态下的桨距角和角速度的数据代入输出向量
Figure BDA0003072478400000131
中，确定工作状态下时滞风电机组桨距子系统的输出向量y(k)，k的取值范围为1至N，k为整数。
设置在k∈[60,115]时时滞风电机组桨距子系统发生参数故障。
图2展示了采用本申请提出的基于滤波的工业时滞系统故障检测方法对时滞风电机组桨距子系统进行故障估计所得到的估计输出的上下界和时滞风电机组桨距子系统的真实输出之间的关系。由图2可以看出，时滞风电机组桨距子系统的真实输出大约在k＝60时超出基于最优全对称多胞体卡尔曼滤波算法估计出的上下界，并且在k＝120左右时回到基于最优全对称多胞体卡尔曼滤波算法估计出的上下界之内，表明时滞风电机组桨距子系统大约在k＝60时发生故障，并且在k＝120左右故障结束，回到正常的工作状态。
图3展示了采用本申请提出的基于滤波的工业时滞系统故障检测方法对时滞风电机组桨距子系统进行故障的故障检测结果。图3所示的仿真曲线显示在k∈[60,120]时时滞风电机组桨距子系统故障指示信号为1，表明时滞风电机组桨距子系统此时发生了故障，故障检测时间与实际故障发生时间接近，表明本申请提出的一种基于滤波的工业时滞系统故障检测方法能够有效对时滞系统进行故障检测的处理，正确检测出时滞系统的故障。
本发明实施例中的部分步骤，可以利用软件实现，相应的软件程序可以存储在可读取的存储介质中，如光盘或硬盘等。
以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种基于滤波的工业时滞系统故障检测方法，其特征在于，所述方法包括：
步骤101，根据工业时滞系统的工作原理，建立工业时滞系统的离散线性模型；
步骤102，获取工业时滞系统在工作状态下的输出向量；
步骤103，根据所建立的工业时滞系统的离散线性模型和工业时滞系统在工作状态下的输出向量，设计针对工业时滞系统的最优全对称多胞体卡尔曼滤波器；
步骤104，根据所设计的针对工业时滞系统的最优全对称多胞体卡尔曼滤波器，估计工业时滞系统的状态向量对应的全对称多胞体；
步骤105，根据状态向量对应的全对称多胞体，计算工业时滞系统的输出向量对应的全对称多胞体的上下界；
步骤106，根据输出向量对应的上下界，确定工业时滞系统的故障情况；
所述步骤101建立的建立工业时滞系统的离散线性模型为：
Figure FDA0003395851570000011
其中x(k)，u(k)和y(k)分别表示工业时滞系统在k时刻的状态向量，输入向量和输出向量，k为离散时间；
w(k)和v(k)分别为过程噪声和测量噪声，并且均有界，分别为w(k)∈W和v(k)∈V，其中W为过程噪声w(k)对应的全对称多胞体，V为测量噪声v(k)对应的全对称多胞体；
A，B，C，D，F分别为对应的参数矩阵，Ah为时滞项参数矩阵，h为时滞时间常数；
所述步骤103包括：
当0≤k≤h时，工业时滞系统的系统时滞对状态估计结果没有影响，按式(2)～(8)设计最优全对称多胞体卡尔曼滤波器：
Figure FDA0003395851570000012
Figure FDA0003395851570000013
Figure FDA0003395851570000014
Figure FDA0003395851570000015
K(k)＝R(k)S-1(k) (6)
L(k)＝AK(k) (7)
Figure FDA0003395851570000016
其中，
Figure FDA0003395851570000017
Figure FDA0003395851570000018
为
Figure FDA0003395851570000019
的降维矩阵，
Figure FDA00033958515700000110
为k时刻状态向量对应的全对称多胞体的生成矩阵，L(k)为k时刻的最优增益矩阵，Gv为测量噪声v(k)对应的全对称多胞体V的生成矩阵，Qv、
Figure FDA0003395851570000021
S(k)、R(k)、K(k)为中间计算变量，“^”表示估计值符号；
当k＞h时，工业时滞系统的系统时滞对状态估计结果存在影响，按式(9)～(20)设计最优全对称多胞体卡尔曼滤波器：
Figure FDA0003395851570000022
Figure FDA0003395851570000023
Figure FDA0003395851570000024
Figure FDA0003395851570000025
Figure FDA0003395851570000026
Figure FDA0003395851570000027
Figure FDA0003395851570000028
Figure FDA0003395851570000029
K1(k)＝R1(k)S-1(k) (17)
K2(k)＝R2(k)S-1(k) (18)
L(k)＝AK1(k)+AhK2(k) (19)
Figure FDA00033958515700000210
其中“∑”表示累加符号，“Π”表示累乘符号，Gw为w(k)对应的全对称多胞体W的生成矩阵，
Figure FDA00033958515700000211
Figure FDA00033958515700000212
为
Figure FDA00033958515700000213
的降维矩阵，
Figure FDA00033958515700000214
为k-h时刻状态向量对应的全对称多胞体的生成矩阵，
Figure FDA0003395851570000031
Figure FDA0003395851570000032
为
Figure FDA0003395851570000033
的降维矩阵，
Figure FDA0003395851570000034
为k-h-i时刻状态向量对应的全对称多胞体的生成矩阵，L(k-i)为k-i时刻的最优增益矩阵，
Figure FDA0003395851570000035
R1(k)、R2(k)、K1(k)、K2(k)为中间计算变量；
所述步骤104包括：
当0≤k≤h时，状态向量对应的全对称多胞体
Figure FDA0003395851570000036
具体可通过式(21)和(22)计算得到：
Figure FDA0003395851570000037
Figure FDA0003395851570000038
其中
Figure FDA0003395851570000039
Figure FDA00033958515700000310
为k+1时刻状态向量对应的全对称多胞体，
Figure FDA00033958515700000311
和
Figure FDA00033958515700000312
为别为
Figure FDA00033958515700000313
的中心和生成矩阵，
Figure FDA00033958515700000314
和
Figure FDA00033958515700000315
分别为k-h时刻状态向量对应的全对称多胞体的中心和生成矩阵，且
Figure FDA00033958515700000316
当k＞h时，状态向量对应的全对称多胞体
Figure FDA00033958515700000317
具体可通过式(23)～(30)计算得到：
Figure FDA00033958515700000318
Figure FDA00033958515700000319
Figure FDA00033958515700000320
Figure FDA00033958515700000321
Figure FDA00033958515700000322
Figure FDA00033958515700000323
Figure FDA00033958515700000324
Figure FDA0003395851570000041
其中H1(k)、H1i(k)、H2(k)、H2i(k)、H3(k)、H3i(k)为中间计算变量，nx为状态向量x(k)的维度，r为全对称多胞体
Figure FDA0003395851570000042
的生成矩阵的维度，
Figure FDA0003395851570000043
为H1(k)的维度，
Figure FDA0003395851570000044
为H2(k)的维度，
Figure FDA0003395851570000045
为H3(k)的维度。
2.根据权利要求1所述的方法，其特征在于，所述步骤105包括：
根据状态向量对应的全对称多胞形
Figure FDA0003395851570000046
基于测量噪声v(k)对应的全对称多胞体V的生成矩阵Gv，按式(31)计算k时刻输出矩阵对应的全对称多胞体
Figure FDA0003395851570000047
Figure FDA0003395851570000048
其中
Figure FDA0003395851570000049
为k时刻锂电池系统输出矩阵对应的全对称多胞体
Figure FDA00033958515700000410
的中心，
Figure FDA00033958515700000411
为k时刻锂电池系统输出矩阵对应的全对称多胞体
Figure FDA00033958515700000412
的生成矩阵，其中
Figure FDA00033958515700000413
表示闵可夫斯基和，“⊙”表示线性映射；
根据输出矩阵对应的全对称多胞体
Figure FDA00033958515700000414
按式(32)和(33)计算k时刻工业时滞系统输出向量对应的上下界：
Figure FDA00033958515700000415
Figure FDA00033958515700000416
其中
Figure FDA00033958515700000417
表述输出向量对应的上界，
Figure FDA00033958515700000418
表示输出向量对应的下界。
3.根据权利要求2所述的方法，其特征在于，所述步骤106包括：
若步骤105得出的输出向量的上下界满足
Figure FDA0003395851570000051
和
Figure FDA0003395851570000052
则表明工业时滞系统未发生故障，此时故障检测信号f(k)＝0；
否则，表明该工业时滞系统在k时刻发生故障，此时故障检测信号f(k)＝1；
其中
Figure FDA0003395851570000053
和
Figure FDA0003395851570000054
分别为输出矩阵
Figure FDA0003395851570000055
和
Figure FDA0003395851570000056
中的元素。
4.根据权利要求3所述的方法，其特征在于，所述工业时滞系统包括时滞风电机组桨距子系统。
5.根据权利要求4所述的方法，其特征在于，所述方法应用于时滞风电机组桨距子系统时，所述步骤101之前还包括：
根据风电机组桨距子系统的工作原理，建立时滞风电机组桨距子系统的系统模型：
Figure FDA0003395851570000057
其中，β为桨距角，βa为桨速度，βr表示桨距参考值，ωn和ζ为时滞风电机组桨距子系统的自然频率和阻尼系数，“·”表示求导符号。
6.根据权利要求5所述的方法，其特征在于，所述时滞风电机组桨距子系统表示为连续时间状态空间方程为：
Figure FDA0003395851570000058
其中x＝[β,βα]T为时滞风电机组桨距子系统的状态向量，u＝βr为时滞风电机组桨距子系统的输入向量，w和v分别表示时滞风电机组桨距子系统的系统的过程噪声和测量噪声。
7.根据权利要求6所述的方法，其特征在于，针对时滞风电机组桨距子系统，
Figure FDA0003395851570000059
Figure FDA0003395851570000061
为时滞项参数矩阵，h＝6为时滞时间常数。
CN202110550989.3A 2021-05-19 2021-05-19 一种基于滤波的工业时滞系统故障检测方法 Active CN113236506B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110550989.3A CN113236506B (zh) 2021-05-19 2021-05-19 一种基于滤波的工业时滞系统故障检测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110550989.3A CN113236506B (zh) 2021-05-19 2021-05-19 一种基于滤波的工业时滞系统故障检测方法

Publications (2)

Publication Number Publication Date
CN113236506A CN113236506A (zh) 2021-08-10
CN113236506B true CN113236506B (zh) 2022-03-15

Family

ID=77137736

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110550989.3A Active CN113236506B (zh) 2021-05-19 2021-05-19 一种基于滤波的工业时滞系统故障检测方法

Country Status (1)

Country Link
CN (1) CN113236506B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117828864B (zh) * 2023-12-29 2024-07-02 江南大学 一种基于多胞空间滤波和p范数的系统状态估计方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2004003910A (ja) * 2002-06-03 2004-01-08 Renesas Technology Corp 半導体集積回路
CN110058124B (zh) * 2019-04-25 2021-07-13 中国石油大学(华东) 线性离散时滞系统的间歇故障检测方法

Also Published As

Publication number Publication date
CN113236506A (zh) 2021-08-10

Similar Documents

Publication Publication Date Title
CN115950557B (zh) 一种基于压力变送器的温度智能补偿方法
CN105843073B (zh) 一种基于气动力不确定降阶的机翼结构气动弹性稳定性分析方法
Abdallah et al. Stator winding inter‐turn short‐circuit detection in induction motors by parameter identification
CN108667673A (zh) 基于事件触发机制的非线性网络控制系统故障检测方法
CN108398637B (zh) 一种非线性机电系统的故障诊断方法
CN113239132B (zh) 一种电压互感器的超差在线辨识方法
CN110838075A (zh) 电网系统暂态稳定的预测模型的训练及预测方法、装置
CN111313405A (zh) 一种基于多量测断面的中压配电网拓扑辨识方法
CN103577710A (zh) 基于分数阶upf的航空功率变换器故障预测方法
CN111046327A (zh) 适用于低频振荡与次同步振荡辨识的Prony分析方法
CN112305418A (zh) 一种基于混合噪声双重滤波的电机系统故障诊断方法
CN113236506B (zh) 一种基于滤波的工业时滞系统故障检测方法
CN115688288B (zh) 飞行器气动参数辨识方法、装置、计算机设备及存储介质
CN113297798A (zh) 一种基于人工神经网络的机器人外界接触力估计方法
CN114372324B (zh) 旋转机械装备关键零部件服役退化轨迹预测方法及设备
CN109462242B (zh) 基于iir数字滤波和esprit辨识算法的电力系统低频振荡检测方法
EP4339621A1 (en) Mems accelerometer systems
JP2867477B2 (ja) オンライン設備の寿命予測方法
Liu et al. An Accelerated Safety Probability Estimation Method for Control Policies of Autonomous Vehicles in Open Environments
CN107732940B (zh) 一种基于adpss的电力系统稳定器参数优化试验方法
JP2021012605A (ja) 伝達関数の予測方法
CN111609878B (zh) 三自由度直升机系统传感器运行状态监测方法
CN111108738A (zh) 数据处理设备、数据分析设备、数据处理系统和用于处理数据的方法
CN108090846B (zh) 一种电网低频振荡案例库的构建方法及装置
Ananthan et al. Novel system model‐based fault location approach using dynamic search technique

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