一种基于集合论和未知输入观测器的故障检测方法

1.本发明涉及故障检测技术，尤其是涉及基于集合论和未知输入观测器的故障检测技术。


背景技术：

2.随着现代工业技术的不断创新与发展，工业生产系统呈现出自动化、规模化与复杂化的特点，与此同时，系统内部发生故障与风险的可能也随之不断增加。因此，具有高可靠性和实时性的、基于系统模型的故障检测技术，对诸如军工生产、核电设施、化工冶金和航空航天等安全性和可靠性要求较高的工业应用场景均具有十分重要的意义。然而，基于模型的故障检测技术依赖于待检测对象的系统模型，但由于扰动、噪声及不确定性因素的存在，在实际建模过程中无法获得真实系统的精确模型。因此，基于模型的故障检测方案必须能够处理扰动、噪声及建模误差对故障检测所带来的影响，这一技术被称为鲁棒故障检测。
3.在20世纪80年代到21世纪初，多数鲁棒故障检测所采用的方案可统称为主动鲁棒故障检测技术，这一时期的故障检测方案大多通过相应的手段去主动解耦噪声、扰动与不确定性因素对故障检测的影响。在此期间，具有代表性的方案有未知输入观测器、特征结构配置、最优奇偶关系以及频域设计方法。然而，主动鲁棒故障检测技术致力于完全解耦噪声、扰动与不确定性的影响，因此上述方案的使用条件较为苛刻，能够应用的场景十分有限。近年来，随着集合论、概率论和数值优化技术的发展，鲁棒故障检测技术发展出了新的范式，即能够容忍噪声、扰动与不确定性因素干扰的被动故障检测技术，此类技术能够在容忍扰动和噪声的前提下实时得出可靠的故障检测阈值，因此在诸多工业生产环节都有着非常广阔的应用前景。
4.现有的被动故障检测技术主要分为两大类：一类是基于扰动和噪声分布已知这一假设的随机故障检测方案，其中代表性的方法有基于卡尔曼滤波器的故障检测技术，此类方案的优点在于能够得到保守性较低的故障检测结果，而缺点是需要得知扰动和噪声分布的先验信息，这一要求在一些大型系统和复杂系统中较难实现，因此应用范围受限。
5.另一类被动故障检测方案则是基于扰动和噪声分布未知但有界的确定故障检测方案，其中代表性的方法有集理论未知输入观测器、集值观测器等，此类方案的优点在于无需事先得知扰动与噪声分布的先验信息，应用范围更广，缺点在于获得的故障检测结果保守性较高，如集理论未知输入观测器在系统主要受测量噪声影响时保守性较高，而集值观测器在系统主要受随机扰动影响时保守性较高。


技术实现要素：

6.本发明的目的在于解决基于扰动和噪声分布未知但有界的此类故障检测方案保守性较高的技术问题，提供一种基于集合论和未知输入观测器的故障检测方法。
7.为实现上述目的，本发明采用以下技术方案：
8.一种基于集合论和未知输入观测器的故障检测方法，包括以下步骤：
9.s1：建立具有有界随机扰动与测量噪声的待检测系统的状态空间模型；
10.s2：根据待检测系统的状态空间模型构造参数待定的未知输入观测器；
11.s3：根据系统模型和观测器模型建立状态估计误差与残差迭代方程；
12.s4：根据迭代方程导出状态估计误差集与残差集迭代公式；
13.s5：建立最小化残差集尺寸的优化问题，求得最小化残差集尺寸的观测器参数；
14.s6：根据观测器参数计算最小残差集，执行故障检测。
15.在一些发明实施例中，步骤s1中的有界随机扰动与测量噪声采用中心对称多面体集合来刻画。
16.在一些发明实施例中，步骤s5中的优化问题为最小化残差集弗罗贝尼乌斯范数尺寸的优化问题。
17.在一些发明实施例中，步骤s1包括以下步骤：
18.s1-1：设置以下形式的线性离散时不变系统：
[0019][0020]
其中，其中，和均为系统的参数矩阵，其中n
x
，ny，nu，n
ω
和n
η
分别为系统状态、系统输出、系统输入、随即扰动与测量噪声向量的维数，且参数矩阵(a，c)可检测。为系统在第k时刻的输入，为系统在第k时刻的输出，为系统在第k时刻的状态，为系统在第k时刻的所受的随机扰动，为系统在第k时刻的测量噪声。
[0021]
s1-2：系统的初始状态，随机扰动以及测量噪声均认为能量有界，并采用中心对称多面体集合来刻画，分别记为和其中对于任意中心对称多面体集合z＝《p，g》，p为该集合的中心向量，g为该集合的生成矩阵，所述中心对称多面体具有以下两点性质：
[0022]
(1)
[0023]
(2)
[0024]
其中，和均为中心对称多面体集合，m表示具有相应维度的线性变换矩阵，算子为闵可夫斯基和，且集合x与y的闵可夫斯基和定义为
[0025]
在一些发明实施例中，步骤s2的方法包括：
[0026]
利用步骤s1-1中的系统(1)中的输入向量uk与输出向量yk，构造如下形式的未知输入观测器：
[0027]
[0028]
其中，n，t，k和h均为参数待定的观测器参数，zk为未知输入观测器在第k时刻的状态，为未知输入观测器在第k时刻对系统的状态估计，为未知输入观测器在第k时刻对系统的输出估计。
[0029]
在一些发明实施例中，步骤s3包括以下步骤：
[0030]
s3-1：定义观测器对系统的状态的估计误差与残差：
[0031][0032]
其中，ek为观测器对系统的状态估计误差，rk为观测器对系统输出的残差信号；
[0033]
s3-2：令未知输入观测器参数k＝k1 k2，根据步骤s3-1中(3)的定义，未知输入观测器对步骤s1-1中的系统(1)中的状态估计误差，可表示为：
[0034][0035]
s3-3：令未知输入观测器参数满足如下约束：
[0036][0037]
则步骤s3-2中式(4)可进一步化为：
[0038]ek 1
＝(a-hca-k1c)ek (e-hce)ω
k-hfη
k 1-k1fηk.
ꢀꢀꢀ
(6)
[0039]
s3-4：根据步骤s3-1中(3)中残差信号定义，未知输入观测器的残差可表示为：
[0040]rk
＝cek fηk.
ꢀꢀꢀ
(7)
[0041]
在一些发明实施例中，，步骤s4包括以下步骤：
[0042]
s4-1：利用步骤s3-3中公式(6)所示的状态估计误差表达式，可获得状态估计误差集的迭代计算公式：
[0043][0044]
其中，ek为第k步系统(1)的状态估计误差集，并假设初始状态估计误差集e0为中心对称多面体集合。
[0045]
s4-2：利用公式(7)所示的残差表达式，可获得残差集的迭代计算公式：
[0046][0047]
其中，rk为第k步系统(1)的残差集。
[0048]
在一些发明实施例中，将采用以下准则来实现鲁棒故障检测：
[0049][0050]
这一准则如图2a和图2b所示，若检测到残差信号rk位于残差集rk外，如图2b所示，则可确认系统发生故障。反之，若检测到残差信号rk落入残差集rk内，如图2a所示，则认为未检测到系统故障。
[0051]
由于残差集rk为包含随机扰动与测量噪声的残差信号集的rk估计集，因此准则(10)仅为系统发生故障的充分条件，换言之，即使系统仍有可能处于故障状态。因此，应设法优化残差集rk的尺寸以求降低故障检测的保守性。在本发明实施例中，采用弗贝尼乌斯范数来度量残差集rk的尺寸：对于中心对称多面体集合其弗罗贝尼乌
斯范数尺寸定义为其中算子||
·
||f表示弗罗贝尼乌斯范数。
[0052]
在一些发明实施例中，步骤s5包括以下步骤：
[0053]
s5-1：建立如下优化问题：
[0054][0055]
其中，hk和k
1，k
分别为未知输入观测器参数h和k1在k时刻的残差集最优弗罗贝尼乌斯范数尺寸值，观测器的其他参数则由公式(5)所述约束确定；
[0056]
s5-2：根据步骤s4-1中状态估计误差集(8)和步骤s4-2中残差集迭代公式(9)，且考虑到系统的初始状态估计误差集e0，随机扰动集w以及测量噪声集v均为中心对称多面体，易得各时刻状态估计误差集ek和残差集rk均为中心对称多面体，将其记为和则k 1时刻状态估计误差集ek和残差集rk的生成矩阵分别可表示为：
[0057][0058][0059]
其中，
[0060][0061][0062]
s5-3：为解决步骤s5-1中的问题(11)，分别定义j＝gw(gw)
t
以及p＝cc
t
，根据步骤s5-2中的式(13)，ψ(r
k 1
)可进一步表示为：
[0063][0064]
其中，
[0065][0066]
考虑到式(14)中tr(fgv(gv)
tft
)项为常量，因此步骤s5-1中的问题(11)可等价为：
[0067][0068]
由于问题(16)为无约束凸优化问题，因此未知输入观测器在第k时刻的增益参数可通过求解以下方程解析求得：
[0069][0070]
其结果为：
[0071][0072]
其中，
[0073]
其中，当且仅当矩阵非奇异时可得到式(18)所示的最优未知输入观测器参数而为四个半正定矩阵(λλ
t
，γγ
t
和)之和，为奇异矩阵当且仅当这四个半正定矩阵的零空间交集除原点外非空。因此在实际系统中，矩阵往往是非奇异矩阵。
[0074]
在一些发明实施例中，步骤s6包括以下步骤：
[0075]
s6-1：将步骤s5-3中式(18)所示的最优未知输入观测器参数表达式代入步骤s5-3中式(15)中，并消去其中的项，可得：
[0076][0077]
式(19)所示的二次矩阵方程即为黎卡提差分方程，当k
→
∞时，有s
∞
＝s
k 1
＝sk，此时可得如下形式的代数黎卡提方程：
[0078][0079]
s6-2：根据代数黎卡提方程求解的相关理论可知，当系统(a，c)可检测时，存在某个s
*
∈{s
∞
}以及相应的l
*
满足：
[0080][0081]
使得矩阵的特征值在单位圆之内，因此当未知输入观测器的参数为l
*
＝[h
∞k1，∞
]时，此时未知输入观测器是稳定的，且l
*
为未知输入观测器在时域无穷远处的稳态最优参数；
[0082]
s6-3：求解步骤s6-1中的代数黎卡提方程(20)获得s
∞
后，将其代入步骤s6-2中式(21)即可获得l
*
＝[h
∞k1，∞
]，其中h
∞
和k
1，∞
即为观测器参数h和k在时域无穷远处的稳态最优值；
[0083]
s6-4：在给定初始状态估计误差集e0并获得h
∞
和k
1，∞
后，将初始状态估计误差集e0与h
∞
和k
1，∞
代入步骤s4-1中式(8)迭代计算各时刻的状态估计误差集ek，进一步根据步骤s4-2中式(9)即可获得各时刻最小frobenius范数下的残差集rk。
[0084]
在实际使用中，步骤s6-1中式(20)所示的黎卡提方程可通过相应的数值方式求解得出。
[0085]
本发明具有如下有益效果：
[0086]
本发明提供了一种统一了集理论未知输入观测器和集值观测器的故障检测方法框架，能够实现更低的保守性，并且对待检测系统所需满足的条件要求更低，适用范围更广。
附图说明
[0087]
图1是本发明实施例一种利用集合论和未知输入观测器的故障检测方法的流程图；
[0088]
图2a是本发明实施例中表示未检测到故障的示意图；
[0089]
图2b是本发明实施例中表示检测到故障的示意图；
[0090]
图3a是本发明实施例中表示在系统输出第一分量上的故障检测示意图；
[0091]
图3b是本发明实施例中表示在系统输出第二分量上的故障检测示意图；
[0092]
图4是本发明实施例的一种利用集合论和未知输入观测器的系统在k＝1到k＝25时刻残差集可视化示意图；
[0093]
图5是本发明实施例与现有的集理论未知输入观测器和集值观测器的在不同扰动
与噪声水平和弗罗贝尼乌斯范数度量下的保守性的对比示意图；
[0094]
图6是本发明实施例与现有的集理论未知输入观测器和集值观测器的在不同时刻和弗罗贝尼乌斯范数度量下的保守性的对比示意图。
具体实施方式
[0095]
以下对本发明的实施方式做详细说明。应该强调的是，下述说明仅仅是示例性的，而不是为了限制本发明的范围及其应用。
[0096]
此外，术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”的特征可以明示或者隐含地包括一个或者更多该特征。在本发明实施例的描述中，“多个”的含义是两个或两个以上，除非另有明确具体的限定。
[0097]
本发明为克服上述背景技术中基于扰动和噪声分布未知但有界的此类故障检测方案保守性较高的缺陷，提供了一种基于集合论和未知输入观测器的故障检测方法。
[0098]
下面结合附图对本发明的实施方式作进一步的介绍说明。需要强调的是，以下所述的实施方式仅为本发明的一种示例性质的实施方式，且本发明并非仅限于该种实施方式。
[0099]
图1所示即为本发明提供的一种利用集合论和未知输入观测器的故障检测实施例，该实施例包括以下步骤：
[0100]
s1：建立具有有界随机扰动与测量噪声的待检测系统的状态空间模型，且有界随机扰动与测量噪声采用中心对称多面体集合来刻画；
[0101]
s2：根据待检测系统的状态空间模型构造参数待定的未知输入观测器；
[0102]
s3：根据系统模型和观测器模型建立状态估计误差与残差迭代方程；
[0103]
s4：根据迭代方程导出状态估计误差集与残差集迭代公式；
[0104]
s5：建立最小化残差集弗罗贝尼乌斯范数尺寸的优化问题，求得最小化残差集尺寸的观测器参数；
[0105]
s6：根据观测器参数计算最小弗罗贝尼乌斯范数残差集，执行故障检测。
[0106]
其中步骤5是观测器参数设计的核心步骤。现有基于集合论和观测器的故障检测技术的差异部分就在于步骤二的观测器结构与步骤五的观测器参数设计方法。本发明最大的贡献在于提出了一个统一现有方法的框架和更低的保守性结果，因此在步骤流程上与现有技术差异不大。
[0107]
以下内容进一步描述本发明实施例。
[0108]
考虑以下形式的线性离散时不变系统：
[0109][0110]
其中，其中，和均为系统的参数矩阵，其中n
x
，ny，nu，n
ω
和n
η
分别为系统状态、系统输出、系统输入、随即扰动与测量噪声向量的维数，且参数矩阵(a，c)可检测。为系统在第k时
刻的输入，为系统在第k时刻的输出，为系统在第k时刻的状态，为系统在第k时刻的所受的随机扰动，为系统在第k时刻的测量噪声。
[0111]
在本实施例中，系统的初始状态，随机扰动以及测量噪声均认为能量有界，并采用中心对称多面体集合来刻画，分别记为中心对称多面体集合来刻画，分别记为和其中对于任意中心对称多面体集合z＝《p，g》，p为该集合的中心向量，g为该集合的生成矩阵。需要说明的是，中心对称多面体具有以下两点性质：
[0112]
(1)
[0113]
(2)
[0114]
其中，和均为中心对称多面体集合，m表示具有相应维度的线性变换矩阵，算子为闵可夫斯基和，且集合x与y的闵可夫斯基和定义为
[0115]
进一步地，利用系统(1)中的输入向量uk与输出向量yk，构造如下形式的未知输入观测器：
[0116][0117]
其中，n，t，k和h均为参数待定的观测器参数，zk为未知输入观测器在第k时刻的状态，为未知输入观测器在第k时刻对系统的状态估计，为未知输入观测器在第k时刻对系统的输出估计。
[0118]
定义观测器对系统的状态的估计误差与残差：
[0119][0120]
其中，ek为观测器对系统的状态估计误差，rk为观测器对系统输出的残差信号。
[0121]
进一步地，令未知输入观测器参数k＝k1 k2，根据(3)中定义，未知输入观测器对系统(1)中的状态估计误差，可表示为：
[0122][0123]
此时，令未知输入观测器参数满足如下约束：
[0124][0125]
则式(4)中所述状态估计误差公式，可进一步化为：
[0126]ek 1
＝(a-hca-k1c)ek (e-hce)ω
k-hfη
k 1-k1fηk.
ꢀꢀꢀ
(6)
[0127]
根据(3)中残差信号定义，未知输入观测器的残差可表示为：
[0128]rk
＝cek fηk.
ꢀꢀꢀ
(7)
[0129]
利用公式(6)所示的状态估计误差表达式，可获得状态估计误差集的迭代计算公
式：
[0130][0131]
其中，ek为第k步系统(1)的状态估计误差集，并假设初始状态估计误差集e0为中心对称多面体集合。进一步地，利用公式(7)所示的残差表达式，可获得残差集的迭代计算公式：
[0132][0133]
其中，rk为第k步系统(1)的残差集。
[0134]
在本发明实施例中，将采用以下准则来实现鲁棒故障检测：
[0135][0136]
这一准则如图2a和图2b所示，若检测到残差信号rk位于残差集rk外，如图2b所示，则可确认系统发生故障。反之，若检测到残差信号rk落入残差集rk内，如图2a所示，则认为未检测到系统故障。
[0137]
由于残差集rk为包含随机扰动与测量噪声的残差信号rk的估计集，因此准则(10)仅为系统发生故障的充分条件，换言之，即使rk∈rk系统仍有可能处于故障状态。因此，应设法优化残差集rk的尺寸以求降低故障检测的保守性。在本发明实施例中，采用弗罗贝尼乌斯范数来度量残差集rk的尺寸：对于中心对称多面体集合其弗罗贝尼乌斯范数尺寸定义为其中算子||
·
||f表示弗罗贝尼乌斯范数。
[0138]
根据上述分析，为降低本发明实施例所述故障检测方法的保守性，可建立如下优化问题：
[0139][0140]
其中，hk和k
1，k
分别为未知输入观测器参数h和k1在k时刻的残差集最优弗罗贝尼乌斯范数尺寸值，观测器的其他参数则由公式(5)所述约束确定。
[0141]
根据状态估计误差集和残差集迭代公式(8)和(9)，且考虑到系统的初始状态估计误差集e0，随机扰动集w以及测量噪声集v均为中心对称多面体，易得各时刻状态估计误差集ek和残差集rk均为中心对称多面体，将其记为和则k 1时刻状态估计误差集ek和残差集rk的生成矩阵分别可表示为：
[0142][0143][0144]
其中：
[0145][0146][0147]
为求解问题(11)，分别定义j＝gw(gw)
t
以及p＝cc
t
。根据式(13)，ψ(r
k 1
)可进一步表示为：
[0148]
[0149]
其中：
[0150][0151]
考虑到式(14)中tr(fgv(gv)
tft
)项为常量，因此问题(11)可等价为：
[0152][0153]
由于问题(16)为无约束凸优化问题，因此未知输入观测器在第k时刻的增益参数可通过求解以下方程解析求得：
[0154][0155]
其结果为：
[0156][0157]
其中，
[0158]
需要说明的是，当且仅当矩阵非奇异时可得到式(18)所示的最优未知输入观测器参数而为四个半正定矩阵(λλ
t
，γγ
t
和)之和，为奇异矩阵当且仅当这四个半正定矩阵的零空间交集除原点外非空。因此在实际系统中，矩阵往往是非奇异矩阵。
[0159]
进一步地，将式(18)所示的最优未知输入观测器参数表达式代入式(15)中，并消去其中的项，可得：
[0160][0161]
而式(19)所示的二次矩阵方程即为黎卡提差分方程。当k
→
∞时，有s
∞
＝s
k 1
＝sk，此时可得如下形式的代数黎卡提方程：
[0162][0163]
根据代数黎卡提方程求解的相关理论可知，当系统(a，c)可检测时，存在某个s
*
∈{s
∞
}以及相应的l
*
满足：
[0164][0165]
使得矩阵的特征值在单位圆之内，因此当未知输入观测器的参数为l
*
＝[h
∞k1，∞
]时，此时未知输入观测器是稳定的，且l
*
为未知输入观测器在时域无穷远处的稳态最优参数。需要说明的是，代数黎卡提方程通常无法解析求解，在实际使用中式(20)所示的黎卡提方程可通过相应的数值方式求解得出。
[0166]
求解步骤s6-1中的代数黎卡提方程(20)获得s
∞
后，将其代入步骤s6-2中式(21)即可获得l
*
＝[h
∞k1，∞
]，其中h
∞
和k
1，∞
即为观测器参数h和k在时域无穷远处的稳态最优值；
[0167]
将初始状态估计误差集e0与h
∞
和k
1，∞
代入步骤s4-1中式(8)迭代计算各时刻的状态估计误差集ek，进一步根据步骤s4-2中式(9)即可获得各时刻最小frobenius范数下的残差集rk。
[0168]
验证例：
[0169]
为验证本发明实施例所述方法的有效性，考虑系统(1)具有如下参数：
[0170][0171][0172]
且随机扰动集w＝《0，0.1i4×4》，测量噪声集v＝《0，0.1i2×2》，未知输入观测器的初始状态估计集e0＝《0，0.01i2×2》。
[0173]
根据本发明实施例所述方法，则未知输入观测器在时域无穷远处的最优参数应设计为：
[0174][0175][0176][0177]
假定系统的初始状态为x0＝[0，0]
t
，输入uk＝[3sin(0.5k)，5sin(0.7k)]
t
，且系统在k＝1到k＝20期间未发生故障，在k＝21时刻，系统执行器发生故障，导致系统的参数矩阵b变为
[0178][0179]
为便于说明，定义rk(1)和rk(2)为残差集rk在第一和第二维分量上的投影，且rk＝[rk(1)，rk(2)]
t
。从图3a和图3b中可以看出，在k＝21时系统发生了故障，未知输入观测器在k＝22时在第二维分量上检测到系统发生了故障，从而证明了本发明实施例的有效性。而图4则为系统在k＝1到k＝25期间残差集rk的可视化示意。
[0180]
比较例：
[0181]
进一步地，为比较本发明所提出故障检测算法和现有的集理论未知输入观测器和集值观测器方法的保守性，考虑系统(1)具有如下参数：
[0182][0183][0184]
且测量噪声集v＝《0，0.1i2×2》，未知输入观测器的初始状态估计集e0＝《0，0.01i2×2》。假定系统的初始状态为x0＝[0，0]
t
，输入uk＝[3sin(0.5k)，5sin(0.7k)]
t
。令随机扰动集的第一分量维度为w1＝《0，αi1×1》，而其余三个分量维度的噪声集为w2＝《0，αi3×3》。则参数α的大小即可反映系统受随机扰动和测量噪声影响水平的相对高低：当参数α较大时，系统主要受随机扰动影响。当参数α较小时，则系统主要受测量噪声影响。
[0185]
为便于说明，将和分别定义为本文所述算法、集理论未知输入观测器和集值观测器三类算法在系统第k步所产生的残差集。则图5即为三种方法在不同扰动与噪声水平和弗罗贝尼乌斯范数度量下的保守性的对比，图中选取的时刻为k＝40(不失一般性，当k＞5时三种方法所产生的残差集均以收敛)，其中α较小时表示系统主要受测量噪声影响，α较大时表示系统主要受随机扰动的影响。可以看出本文所述算法的保守性最低。进
一步地，图6为三类方法在不同时刻和弗罗贝尼乌斯范数度量下的保守性的对比，图中选取的α＝0.2，图6同样验证了本文所述算法在三类算法中保守性最低。
[0186]
以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的技术人员来说，在不脱离本发明构思的前提下，还可以做出若干等同替代或明显变型，而且性能或用途相同，都应当视为属于本发明的保护范围。