-
摘要: 对于现代复杂控制系统, 微小故障往往很难发现. 在系统过程干扰和测量噪声未知但有界的前提下, 提出了一种新的基于状态集员估计的主动故障检测方法. 首先设计全对称多胞形卡尔曼滤波器对系统状态进行估计, 并利用全对称多胞形对受未知干扰影响的状态集合进行描述, 然后设计辅助输入信号使得加入辅助输入信号后正常模型的状态集合与故障模型的状态集合交集为空, 从而实现主动故障检测. 为了使得所设计的辅助输入信号对原系统影响最小, 需要求得最小的辅助输入信号, 本文将最优化问题转化为混合整数二次规划问题进行求解. 最后, 与基于输出集合的辅助输入信号设计方法对比, 仿真验证本文所提出的基于状态集合的主动故障检测方法由于未受下一时刻测量噪声的影响, 所求得的辅助输入信号更小, 保守性更低.Abstract: For modern complex systems, incipient faults is usually difficult to be detected. Under the assumption that system process disturbance and measurement noise are unknown but bounded, this paper proposes a novel active fault detection method based on state set-membership estimation. First, a zonotopic Kalman filter is designed to estimation systems states, and the state sets affected by the unknown inputs are described by zonotopes. Then, an auxiliary input signal is designed such that the state sets of the normal model are separated from the ones of faulty models, as a result, incipient faults are detected successfully. In order to decrease the effect of the auxiliary input to the practical systems, the minimum auxiliary input signal is required. In this paper, the optimization problem is transformed into a mixed integer quadratic programming problem. Compared with the output sets based auxiliary input signal design method, the proposed technique can achieve a smaller auxiliary input signal because of states set are not affect by the measurement noise at the next time instant, and it has less conservatism.
-
Key words:
- Active fault detection /
- state set-membership estimation /
- auxiliary input signal /
- zonotope /
- incipient fault
-
系统的操作安全以及产品的质量对于工业过程控制至关重要, 而一旦系统发生故障不仅会影响产品质量, 甚至会给整个系统带来安全隐患, 因此故障检测及诊断具有重要的研究意义与研究价值[1]. 通常, 系统根据是否发生故障及故障的大小可分为三个阶段: 第1阶段为正常阶段, 此时系统虽然受到未知干扰和测量噪声的影响, 但是被监控系统的变量处在可接受的范围内. 当系统的被监控变量的幅值或者特征发生微小变化, 使得系统的某些量超出可接受的阈值, 但是仍小于某一较大的阈值时, 此时为第2阶段. 这些微小的变化一般很难被观测到, 此阶段为微小故障或者潜伏故障阶段, 该阶段的变化不会大程度地影响系统. 但是微小故障会随着时间的变化逐渐变大, 当超过较大的这个阈值时, 此时系统处在第3阶段, 即显著故障阶段. 该阶段被监控变量的幅值发生大幅度变化, 很容易被监控人员发现[2]. 系统处在第3阶段时, 由于系统的故障特征十分明显, 往往会存在较为重大的安全隐患. 因此, 针对此阶段的故障检测与诊断引起了高度重视. 在此阶段通过系统的测量数据生成残差, 而无需设计辅助输入信号激励系统来判断故障发生与否, 此类故障诊断方法通常称为被动故障诊断. 被动故障诊断主要分为三类: 基于模型的方法[3-5]、基于知识的方法[6-8]以及基于数据驱动[9-11]的方法.
当系统处在第2阶段, 即系统发生微小故障时, 由于其对系统的影响较小, 往往未造成实质性伤害, 此时, 对微小故障进行及时有效地监控, 将对预防和减少灾难的发生具有重要作用[12]. 但同时由于微小故障发生在初始阶段, 幅值较小, 受未知干扰信号的影响, 采用被动故障检测及诊断方法往往不能准确地检测出来. 相较于被动故障诊断, 主动故障诊断方法通过生成辅助输入信号对原系统加以激励, 从而使得微小故障的影响得以显现, 近年来受到越来越多研究者的关注[13-14]. 主动故障诊断方法根据未知干扰信号的处理方式可以分为确定性方法和概率性方法[15]. 确定性方法是将系统的干扰和噪声等不确定的因素假设为具有已知上下界或者已知内部结构分布的信号, 如控制器和检测器参数设计[16-17]、设计辅助输入信号[18]等方法. 概率性方法是将系统不确定性假设为具有已知概率密度分布的变量, 如结合贝叶斯方法的主动故障诊断[19]等方法.
由于基于辅助输入信号设计的主动故障诊断原理较为简单, 且不会改变系统的结构, 所以受到研究者关注. 传统基于辅助输入信号设计的主动故障诊断方法一般假设系统的干扰和测量噪声服从某种概率分布, 而在实际工业中的干扰和噪声分布的先验信息往往很难获得. 为此, 研究者假设干扰和噪声内部结构未知但是上下界已知, 提出了集员估计理论[20]. 由于该假设符合实际工业过程干扰和噪声的特性, 得到研究者的青睐. 集员估计理论期望得到与系统模型、测量输出以及未知干扰相包容的状态或参数的可行集, 通常利用几何体, 如椭球、平行多面体、区间、全对称多胞形(Zonotope)等来近似该可行集[21]. 文献[22]针对于多模型系统, 利用不同模型输出所在椭球交集为空, 设计最小辅助输入信号, 对不同模型类型进行检测分离. 文献[23]针对线性时不变系统, 利用0~$N$时刻输出所在的椭球, 实现了离线故障检测及辅助输入信号设计. 相比较于椭球, 全对称多胞形的结构以及运算方式简单, 保守性较小等优点, 近几年得到大量研究. 针对线性加性故障系统, 文献[24] 通过引入全对称多胞形的思想, 设计并比较了开环输入信号和闭环输入信号的大小. 文献[25]利用0~$N$时刻正常和故障模型输出所在的全对称多胞形的交集为空, 进行故障检测, 同时得到相应辅助输入信号激励系统.
主动故障检测主要是利用设计得到的辅助输入信号激励系统, 实现微小故障检测. 因此在保证该辅助输入信号能够激励故障的前提下, 同时也应该对系统的影响较小. 目前, 基于集员估计的主动故障检测的研究主要集中在通过分析系统的输出所在的集合, 来进行故障检测和辅助输入信号设计[26]. 但是由于系统输出的测量值所在的集合中不仅包含了系统的过程干扰也包含了测量噪声等不确定因素, 其保守性较大. 而系统的状态所在的集合比输出所在的集合少了下一时刻的测量噪声, 因此利用状态所在集合交集为空设计辅助输入信号可降低一部分保守性, 使得设计得到的辅助输入信号较小. 为此, 本文主要针对微小故障, 研究基于状态集员估计的辅助输入信号设计方法. 首先, 利用全对称多胞形卡尔曼滤波对系统状态进行估计. 然后利用该观测器观测得到系统状态所在的集合, 在保证辅助输入信号最小的情况下, 通过使正常和故障模型的状态所在的全对称多胞形的交集为空, 得到优化问题. 将该优化问题转化为混合整数二次规划(Mixed integer quadratic programming, MIQP), 通过求解该MIQP问题, 得到最优的辅助输入信号.
1. 问题描述及相关概念
1.1 问题描述
考虑如下线性离散不确定系统
$$ \begin{align} \begin{cases} {\pmb x}_{k+1}^{[i]}=A^{[i]}{\pmb x}_{k}^{[i]}+{{B}^{[i]}}{{\pmb u}_{k}}+{{E}^{[i]}}{{\pmb w}_{k}^{[i]}} \\[1mm] {\pmb y}_{k}^{[i]}={{C}^{[i]}}{\pmb x}_{k}^{[i]}+{{F}^{[i]}}{{\pmb v}_{k}^{[i]}} \end{cases} \end{align} $$ (1) 其中, $i \in \{0, 1, \cdots, q\}$, ${\pmb x}_k\in{{\bf{R}}^{n_x}}$, ${\pmb u}_k\in{{\bf{R}}}^{n_u}$, ${\pmb y}_k\in$ ${{\bf{R}}^{n_y}}$ 分别为系统在$k$ 时刻的状态, 输入及测量输出向量. $q$ 代表发生故障的个数. 当$i=0$ 时, 代表系统正常, 当$i$ 为其他值时, 则代表系统发生故障. $A^{[i]}$, $B^{[i]}$, $C^{[i]}$, $E^{[i]}$, $F^{[i]}$ 分别为适当维度的参数矩阵. ${\pmb w}_k\in{{\bf{R}}^{n_w}}$ 和${\pmb v}_k\in{{\bf{R}}^{n_v}}$ 分别为系统的过程干扰及测量噪声. 不失一般性, 本文假设系统的初始状态, 过程干扰和测量噪声分别满足: ${\pmb x}_0 \in \mathcal{X}_0$, ${\pmb w}_k\in$ $\mathcal{W}$, ${\pmb v}_k \in \mathcal{V}$, 其中, $\mathcal{X}_0$, $\mathcal{W}$, $\mathcal{V}$ 分别为初始状态、过程干扰以及测量噪声所在的全对称多胞形集合.
通常, 微小故障很难被观测人员察觉, 主动故障检测利用系统的输出数据产生决策及辅助输入信号, 将该辅助输入信号注入给系统, 提高检测质量. 集员估计理论是主动故障检测常用的一种方法, 将系统的干扰及测量噪声等不确定的因素用集合的形式表示成一种几何形状, 如区间、椭球、多面体、全对称多胞形等, 然后设计最优的辅助输入信号, 使得正常模型和故障模型的输出所在的集合交集为空, 即
$$ \begin{align} &{{\mathcal{Y}}^{[0]}} \cap {{\mathcal{Y}}^{[i]}}=\emptyset, \quad ~i \in \{1, \cdots, q\} \end{align} $$ (2) 其中, $\mathcal{Y}^{[0]}$ 代表正常模型输出所在的集合, ${\mathcal{Y}}^{[i]}$ 为故障模型输出所在的集合.
主动故障检测的动态过程如图 1所示. 其中, 图 1 (a)是系统的辅助输入信号为零时输出的集合, 当故障为微小故障时两个集合相交, 当实际输出属于相交部分时无法确认系统是正常还是故障. 基于输出集合的主动故障检测的原理是通过设计辅助输入信号, 使得正常输出集合与故障模型输出集合相分离, 如图 1 (b)所示.
图 1 在正常模型和故障模型下, 加入辅助信号与未加辅助信号结果的对比Fig. 1 Comparison of the results of adding auxiliary signals and unassisted signals under normal and fault models注意到已有的方法采用的输出集合交集为空的方法, 由式(1)可得输出集合的表示
$$ \begin{align} \mathcal{Y}=C\mathcal{X}+F\mathcal{V} \end{align} $$ (3) 其中, $\mathcal{Y}$ 代表输出所在的集合, $\mathcal{X}$ 为状态所在的集合, $\mathcal{V}$ 为测量噪声所在的集合.
由于$C\mathcal{X}$ 为线性变化, 并不影响两个交集的相交与否, 如图 2 (b)所示. 但由于$F\mathcal{V}$ 的引入, 使得原本相切的两个集合相交, 相交面积的大小与引入的测量噪声$\mathcal{V}$ 的大小有关. 若针对图 2 (c)中利用输出集合设计辅助输入信号使得输出集合交集为空, 则引入的信号必将比利用状态集合的交集为空的方法大, 进而对原系统的影响大.
图 2 在正常模型和故障模型下, 加入辅助信号对状态和输出的影响Fig. 2 The effect of adding an auxiliary signal on state and output under normal and fault models综上分析, 利用状态交集为空设计辅助输入信号的保守性较低. 因此, 本文主要针对系统(1) 的状态所在的集合交集为空设计最优辅助输入信号, 用该最优辅助输入信号激励系统, 提高故障检测质量.
1.2 基本概念
目前主动故障检测方法中常用的集合有椭球、全对称多胞形等. 相较于其他集合, 全对称多胞形的形式简单、计算灵活性高, 可以描述任意凸面体, 并能够适当的控制围包效应[27]等. 因此, 本文主要研究利用全对称多胞形对状态集合进行描述, 相关定义与基本性质描述如下.
定义1[25]. $r$ 阶全对称多胞形$\mathcal{Z}$ 的表达式定义为
$$ \begin{align} \mathcal{Z}=\pmb c+\sum\limits_{i=1}^{r}{{{\alpha }_{i}}{{\pmb g}_{i}}}= \left\langle \pmb c, G \right\rangle =\pmb c\oplus G{{\pmb {B}}^{r}} \end{align} $$ (4) 其中, $-1\le {{\alpha }_{i}}\le 1$, 超立方体${{\pmb {B}}^{r}}={{\left[-1, +1 \right]}^{r}}$, $\pmb c$ $\in$ ${\bf{R}}^n$ 为$\mathcal{Z}$ 的中心, 向量${\pmb g_1, \pmb g_2, \cdots, \pmb g_r}$ 为$\mathcal{Z}$ 的生成器, $G={\{\pmb g_1, \pmb g_2, \cdots, \pmb g_r\}}\in {\bf{R}}^{n\times r}$ 为$\mathcal{Z}$ 的生成器矩阵.
全对称多胞形满足如下基本计算规则:
$$ \left\langle {{\pmb c}_{1}}, {{G}_{1}} \right\rangle \oplus \left\langle {{\pmb c}_{2}}, {{G}_{2}} \right\rangle =\left\langle {{\pmb c}_{1}}+{{\pmb c}_{2}}, [{{G}_{1}}\text{ }{{G}_{2}}] \right\rangle $$ (5) $$ M\otimes \left\langle \pmb c, G \right\rangle =\left\langle M{\pmb c}, MG \right\rangle $$ (6) $$ -\mathcal{Z}=\left\langle -\pmb c, G \right\rangle $$ (7) Frobenius半径可以作为衡量全对称多胞形大小的指标, 一般由生成矩阵$G$的Frobenius范数表示:
$$ \begin{align} & {{J}_{\rm F}}=\sqrt{{\rm tr}({{G}^{\rm T}}G)}=\sqrt{{\rm tr}(G{{G}^{\rm T}})} \end{align} $$ (8) 由式(5)可知, 随着时间的推移, 全对称多胞形的阶数会不断增加, 最终无法计算, 为此, 需要对全对称多胞形的生成矩阵的维数加以限制.
引理1[28].全对称多胞形可以由一个最小的盒子包含
$$ \begin{align} \mathcal{Z}\subset \square \mathcal{Z}=\pmb c\oplus rs(G){{\pmb {B}}^{r}} \end{align} $$ (9) 其中, $\square \mathcal{Z} $代表包含集合$\mathcal{Z}$ 的最小盒子, $rs{(G)}$ 为对角矩阵, 满足
$$ \begin{align} %rs{{(G)}_{i, j}}=\sum\nolimits_{j=1}^{r}{\left| {{G}_{i, j}} \right|} rs{{(G)}}=\left[ \begin{matrix} \sum\limits_{j=1}^{r}{\left| {{G}_{1, j}} \right|} & 0 & \cdots & 0 \\ 0 & \sum\limits_{j=1}^{r}{\left| {{G}_{2, j}} \right|} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \sum\limits_{j=1}^{r}{\left| {{G}_{n, j}} \right|} \\ \end{matrix} \right] \end{align} $$ (10) 引理2[29].考虑$r$ 维的全对称多胞形$\mathcal{Z}= \pmb c $ $\oplus$ $G\pmb{B}^{r} \subset {\bf{R}}^n $ 以及整数$s$, 且满足$n\leq s \leq r$, 定义$\bar{G}$ 为将$G$ 的列向量按欧氏范数降序排列所得矩阵, 则全对称多胞形$ \mathcal{Z} $ 可以包含一个$ s $ 维的全对称多胞形, 即
$$ \begin{align} \mathcal{Z}\subseteq \pmb c\oplus \begin{bmatrix}\bar{G}_1 &Q\end{bmatrix} B^s \end{align} $$ (11) 其中, $\bar{G}_1$ 由$\bar{G}$ 的前$ (s-n) $ 个列向量组成, 定义$\bar{G}_2$ 为$\bar{G}$ 中除去$\bar{G}_1 $ 的部分, $Q \pmb{B}^n $ 为按引理1求得的包含$\bar{G}_2\pmb{B}^{r-s-n} $ 的区间包络, 即
$$ \begin{align} Q_{tt}=\sum\limits_{j=1}^{r} \vert \bar{G}_{2}\vert_{tj}, \quad t=1, \cdots, n \end{align} $$ (12) 引理3[25].已知全对称多胞形$\mathcal{Z}=\langle {{{\pmb a}}_{z}}+{{{\pmb b}}_{z}}$, ${{{G}}_{z}} \rangle $ 和$\mathcal{Y}= \langle {{\pmb a}_{y}}+{ {\pmb b}_{y}}, {{G}_{y}} \rangle $, 则$\mathcal{Z}\cap \mathcal{Y}=\emptyset$ 的充要条件为
$$ \begin{align} {{\pmb a}_{y}}-{ {\pmb a}_{z}}\notin \left\langle {{\pmb b}_{z}}, {{G}_{z}} \right\rangle \oplus \left\langle { -{\pmb b}_{y}}, { {G}_{y}} \right\rangle \end{align} $$ (13) 2. 全对称多胞形卡尔曼滤波器设计
本文主要研究基于状态集合交集为空设计辅助输入信号的方法. 由于状态量无法实时获取, 本节设计一个全对称多胞形卡尔曼滤波器对状态集合进行估计, 并利用全对称多胞形对受干扰影响的状态集合进行描述.
假设系统(1)的初始状态, 过程干扰和测量噪声分别满足:
$$ \begin{align} & {{\pmb x}_{0}^{[i]}}\in {\mathcal{X}_{0}^{[i]}}=\left\langle {{\pmb c}_{0}^{[i]}}, {{G}_{0}^{[i]}} \right\rangle \nonumber\\[0.5mm] & {{\pmb w}_{k}^{[i]}}\in \mathcal {W}^{[i]}=\left\langle 0, {{G}_{w}^{[i]}} \right\rangle \nonumber\\[0.5mm] & {{\pmb v}_{k}^{[i]}}\in \mathcal {V}^{[i]}=\left\langle 0, {{G}_{v}^{[i]}} \right\rangle \end{align} $$ (14) 针对该系统, 设计全对称多胞形卡尔曼滤波器形式如下:
$$ \begin{align} \hat{\pmb x}_{k+1}^{[i]}=&\ {{A}^{[i]}}\hat{\pmb x}_{k}^{[i]}\, + {{B}^{[i]}}{{\pmb u}_{k}}+{{E}^{[i]}}\pmb w_{k}^{[i]}\, +\nonumber\\ & \ L_{k}^{[0]}(\pmb y_{k}^{[i]}-{{C}^{[i]}}\hat{\pmb x}_{k}^{[i]}-{{F}^{[i]}}\pmb v_{k}^{[i]}) \end{align} $$ (15) 其中, $L_{k}^{[0]}$为滤波器的增益矩阵, $\hat{\pmb x}_{k}^{[i]}$ 为$k$ 时刻状态的估计值.
定理1. 设初始时刻的状态满足$\pmb x_{0}^{[i]}\in \mathcal{X}_{0}^{[i]}$, 已知$k$ 时刻估计得到的状态为$\hat{\pmb x}_{k}^{[i]}\in \hat{\mathcal{X}}_{k}^{[i]}=\langle \pmb c_{k}^{[i]}, G_{k}^{[i]}\rangle$, 则$k+1$ 时刻的状态为
$$ \begin{align} \hat{x}_{k\text{+}1}^{[i]}\in \hat{\mathcal{X}}_{k\text{+}1}^{[i]}=\left\langle \pmb c_{k\text{+}1}^{[i]}, G_{k+1}^{[i]} \right\rangle \end{align} $$ (16) 其中,
$$ \begin{align} &\pmb c_{k+1}^{[i]}=({{A}^{[i]}}-{L_{k}^{[0]}}{{C}^{[i]}})\pmb c_{k}^{[i]}+{{B}^{[i]}}\pmb u_{k}^{[i]}+{L_{k}^{[0]}}\pmb y_{k}^{[i]} \\[2mm] &G_{k+1}^{[i]}\!=\! [ \begin{matrix} ({{A}^{[i]}}-{L_{k}^{[0]}}{{C}^{[i]}}){\bar G}_{k}^{[i]} &\! {{E}^{[i]}}G_{w}^{[i]} & \!-{L_{k}^{[0]}}{{F}^{[i]}}G_v^{[i]} \\ \end{matrix} ] \end{align} $$ ${\pmb c_{k+1}^{[i]}}\in {\bf R}^n$, ${G_{k+1}^{[i]}}\in {\bf R}^{n\times r}$, ${\bar G}_{k}^{[i]}$表示将$k$时刻状态所在全对称多胞形降阶后的生成矩阵.
证明.由式(15)可得:
$$ \begin{align} \hat{\pmb x}_{k+1}^{[i]} \in \hat{\mathcal{X}}_{k\text{+}1}^{[i]} =&\left[ {{A}^{[i]}}\odot \left\langle \pmb c_{k}^{[i]}, G_{k}^{[i]} \right\rangle \right]\oplus {{B}^{[i]}}{{\pmb u}_{k}}\oplus\nonumber\\ & \left[ {{E}^{[i]}}\odot \left\langle 0, {{G}_{w}^{[i]}} \right\rangle \right]\oplus {L_{k}^{[0]}}\pmb y_{k}^{[i]} \oplus\nonumber\\ & \left[ (-{L_{k}^{[0]}}{{C}^{[i]}})\odot \left\langle \pmb c_{k}^{[i]}, G_{k}^{[i]} \right\rangle \right] \oplus\nonumber\\ & \left[ (-{L_{k}^{[0]}}{{F}^{[i]}})\odot \left\langle 0, {{G}_{v}^{[i]}} \right\rangle \right] \end{align} $$ (17) 根据全对称多胞形的性质以及引理2, 化简式(17)可得$k+1$时刻估计得到的状态所在的全对称多胞形的中心$ \pmb c_{k+1}^{[i]}$ 及生成矩阵$G_{k+1}^{[i]}$.
定理2. 依据卡尔曼滤波器设计方法, 针对状态全对称多胞形(16), 设计的增益矩阵为式(18)~(23), 则观测器(15)可以实时估计系统(1)的状态所在的集合.
$$ \begin{align} & L_{k}^{[0]*}\text{ } = {{A}^{[0]}}K_{k}^{[0]*}\text{ } \end{align} $$ (18) $$ \begin{align} K{{_{k}^{[0]*}}} = R_{k}^{[0]}(S_{k}^{[0]})^{-1}\text{ } \end{align} $$ (19) $$ \begin{align} R_{k}^{[0]} = \bar{P}_{k}^{[0]}({C}^{[0]})^{\rm T}\text{ } \end{align} $$ (20) $$ \begin{align} S_{k}^{[0]} = {{C}^{[0]}}\bar{P}_{k}^{[0]}({C}^{[0]})^{\rm T}+Q_{v}^{[0]}\text{ } \end{align} $$ (21) $$ \begin{align} \bar{P}_{k}^{[0]} = \bar{G}_{k}^{[0]}(\bar{G}_{k}^{[0]})^{\rm T}\text{ } \end{align} $$ (22) $$ \begin{align} Q_{v}^{[0]} = {{F}^{[0]}}({F}^{[0]})^{\rm T}\text{ } \end{align} $$ (23) 证明.定义$J_{{{F}_{k}}}^{2}={\rm tr}(G_{k}^{[0]}(G_{k}^{[0]})^{\rm T})$ 为评判全对称多胞形大小的指标, 当$\frac{\partial J_{{{F}_{k}}}^{2}}{\partial {L_{k}^{[0]}}}=0$时得到最优的观测器增益${L_{k}^{[0]}}=L_{k}^{[0]*}$. 具体证明过程详见文献[30], 本文不再赘述.
3. 最优辅助输入信号设计
本节利用估计的状态集合来设计辅助输入信号, 通过引入辅助输入信号去激励实际系统, 从而实现故障检测. 由于辅助输入信号的加入很可能会影响系统的实际运行状态, 需要设计的辅助输入信号越小越好. 本节将辅助输入信号的设计转化为求解混合整数二次规划问题, 得到最优的辅助输入信号.
为了简化后文计算过程的表达式, 定义
$$ \tilde{\pmb c}_{k+1}^{[i]}= ({{A}^{[i]}}-{L_{k}^{[0]}}{{C}^{[i]}}) \pmb c_{k}^{[i]}+{L_{k}^{[0]}}\pmb y_{k}^{[i]} $$ (24) $$ \begin{align} {\mathcal M}_{k+1}^{[i]}=&\ {{A}^{[i]}}\hat{\mathcal X}_{k}^{[i]} \oplus {{E}^{[i]}}{\mathcal W}^{[i]}\oplus {L_{k}^{[0]}}\pmb y_{k}^{[i]}\, \oplus \nonumber\\[1mm] &\ (-{L_{k}^{[0]}}{{C}^{[i]}}\hat{\mathcal X}_{k}^{[i]}) \oplus (-{L_{k}^{[0]}}{{F}^{[i]}}{\mathcal V}^{[i]})= \nonumber\\[1mm] &\, \left\langle \tilde{\pmb c}_{k+1}^{[i]}, G_{k+1}^{[i]} \right\rangle \end{align} $$ (25) 定理3. 针对系统(1), 设计辅助输入信号${\pmb u}_k$, 满足
$$ \begin{align} & \vartriangle ^{[0i]}{{\pmb u}_{k}}\notin \mathcal{Z}_{m}(0, i) \end{align} $$ (26) 则在辅助输入信号${\pmb u}_k$影响下的正常模型与故障模型状态集合的交集为空, 即
$$ \begin{align} &\hat{\mathcal{X}}_{k+1}^{[0]}\cap \hat{\mathcal{X}}_{k+1}^{[i]}=\emptyset, \quad i\in\{1, 2, \cdots, q\} \end{align} $$ (27) 其中, $\vartriangle ^{[0i]}={{B}^{[i]}}-{{B}^{[0]}}$, $\mathcal{Z}_m(0, i)={{\mathcal{M}}^{[0]}}\oplus (-{{\mathcal{M}}^{[i]}})$.
证明.由定理1以及式(27)整理得:
$$ \begin{align} \hat{\mathcal X}_{k+1}^{[i]}=&\ {{A}^{[i]}}\hat{\mathcal X}_{k}^{[i]} \oplus {{B}^{[i]}}{{\pmb u}_{k}^{[i]}}\oplus {{E}^{[i]}}{\mathcal W}^{[i]}\, \oplus\nonumber\\[1mm] &\ {L_{k}^{[0]}}\pmb y_{k}^{[i]} \oplus (-{L_{k}^{[0]}}{{C}^{[i]}} \hat{\mathcal X}_{k}^{[i]})\, \oplus \nonumber\\[1mm] &\ (-{L_{k}^{[0]}}{{F}^{[i]}}{\mathcal V}^{[i]})=\nonumber\\[1mm] &\ \left\langle \tilde{\pmb c}_{k+1}^{[i]}+B^{[i]}\pmb u_{k}^{[i]}, G_{k+1}^{[i]}\right\rangle \end{align} $$ (28) 为了实现$\hat{\mathcal{X}}_{k+1}^{[0]}\cap \hat{\mathcal{X}}_{k+1}^{[i]}=\emptyset $, 即
$$ \begin{align} &\left\langle \tilde{\pmb c}_{k+1}^{[0]}+B^{[0]}\pmb u_{k}^{[0]}, G_{k+1}^{[0]}\right\rangle \cap\nonumber\\[1mm]& \qquad\left\langle \tilde{\pmb c}_{k+1}^{[i]}+B^{[i]}\pmb u_{k}^{[i]}, G_{k+1}^{[i]}\right\rangle =\emptyset \end{align} $$ (29) 根据引理3, 需要满足
$$ \begin{align} & {{B}^{[i]}}{{\pmb u}_{k}}-{{B}^{[0]}}{{\pmb u}_{k}}\notin \left\langle \tilde{\pmb c}_{k+1}^{[0]}, G_{k+1}^{[0]} \right\rangle \oplus \nonumber\\[1mm] & \qquad\left\langle -\tilde{\pmb c}_{k+1}^{[i]}, G_{k+1}^{[i]} \right\rangle {{\mathcal{M}}^{[0]}}\oplus (-{{\mathcal{M}}^{[i]}}) \end{align} $$ (30) 经整理可得式(27).
为了使得所设计的辅助输入信号对系统的影响最小, 需所设计的辅助输入信号最小, 因此通过求解以下优化问题, 得到辅助输入信号:
$$ \begin{align} &\min\, {\pmb u}_{k}^{\rm T} { R}{{\pmb u}_{k}}\notag\\ &\ {\rm s.t.}\quad {{\vartriangle }^{[0i]}}{{\pmb u}_{k}}\notin \mathcal{Z}_m(0, i), ~~\forall i\in \{1, 2, \cdots , q\} \end{align} $$ (31) 其中, $R$为半正定矩阵. 依据文献[25], 由于优化问题(31)是一个非凸优化问题, 不易得到最优解, 为了得到有效的求解方法, 结合定理3, 本文对优化问题(31)进行重构, 将式(31)转化为混合整数二次规划问题进行求解.
设全对称多胞形${\mathcal{Z}}_m(0, i)=\langle {\pmb c_{m}^{[0i]}}, {G_{m}^{[0i]}} \rangle$, 其中, ${\pmb c_{m}^{[0i]}}$ $\in {\bf R}^n$, $G_{m}^{[0i]}\in {\bf R}^{n\times 2r}$. 当$G_m^{[0i]}$ 行满秩时, 则意味着$\mathcal{Z}_m(0, i)$ 是非空的集合. 当${{\vartriangle }^{[0i]}}{{\pmb u}_{k}}\in \mathcal{Z}_m(0, i)$ 时
$$ \begin{align} {{\vartriangle }^{[0i]}}{{\pmb u}_{k}}={G_{m}^{[0i]}}{{\pmb \xi^{[0i]} }}+{\pmb c_{m}^{[0i]}} \end{align} $$ (32) 其中
$$ \begin{align} {{\left\|\pmb \xi^{[0i]} \right\|}_{\infty }}\le 1 \end{align} $$ (33) 因此, 当${{\vartriangle }^{[0i]}}{{\pmb u}_{k}}\notin \mathcal{Z}_m(0, i)$时
$$ \begin{align} {{\vartriangle }^{[0i]}}{{\pmb u}_{k}}={G_{m}^{[0i]}}{{\pmb \xi^{[0i]} }}+{\pmb c_{m}^{[0i]}} \end{align} $$ (34) 其中
$$ \begin{align} {{\left\|\pmb \xi^{[0i]} \right\|}_{\infty }}\le 1+{{\delta }^{[0i]}} \end{align} $$ (35) 当${{\delta }^{[0i]}}\le 0$ 时, ${{\left\|\pmb \xi^{[0i]} \right\|}_{\infty }}\le 1$, 则不满足${{\vartriangle }^{[0i]}}{{\pmb u}_{k}}$ $\notin$ $\mathcal{Z}_m(0, i)$ 的条件, 只有当${{\delta }^{[0i]}}> 0$ 时, 才能满足该条件.
因此对于每个$\pmb u_k$, 则可以将优化问题(31)重构为线性规划问题
$$ \begin{align} &\min\, {\pmb u}_{k}^{\rm T} { R}{{\pmb u}_{k}}\nonumber\\[2mm] &\rm {s.t.}\quad \begin{cases} {\hat \delta ^{[0i]}} = \mathop {\min }\limits_{{\delta ^{[0i]}},{{\pmb \xi }^{[0i]}}}, {{\delta }^{[0i]}}\\[1mm] {{\vartriangle }^{[0i]}}{{\pmb u}_{k}}={G_{m}^{[0i]}} {{\pmb \xi }^{[0i]}}+{\pmb c_{m}^{[0i]}} \\[1mm] {{\left\| {{\pmb \xi }^{[0i]}} \right\|}_{\infty }}\le 1+{{\delta }^{[0i]}}\\[2mm] {{\hat{\delta }}^{[0i]}}>0 \end{cases} \end{align} $$ (36) 当$G_{m}^{[0i]}$不是行满秩矩阵时, 全对称多胞形$\mathcal{Z}_m{(0, i)}$是空的集合, 则为了使得$\mathcal{Z}_m{(0, i)}$为非空集合, 引入$H_{m}^{[0i]}$矩阵, 使得矩阵$[\begin{matrix} {G_{m}^{[0i]}}~{H_{m}^{[0i]}} \\ \end{matrix}]$行满秩的同时矩阵$H_{m}^{[0i]}$的列数最小. 因此优化问题(36)中的约束条件被替换为
$$ \begin{align} &{{\Delta }^{[0i]}}{{\pmb u}_{k}}={G_{m}^{[0i]}}{{\pmb \xi }^{[0i]}}+{H_{m}^{[0i]}}{{\pmb \gamma }^{[0i]}}+{\pmb c_{m}^{[0i]}}\nonumber\\[1mm] &{{\left\| {{\pmb \xi }^{[0i]}} \right\|}_{\infty }}\le 1+ {{\delta }^{[0i]}}\nonumber\\[1mm] &{{\left\| {{\pmb \gamma }^{[0i]}} \right\|}_{\infty }}\le {{\delta }^{[0i]}} \end{align} $$ (37) 由式(36)中的约束条件${{\hat{\delta }}^{[0i]}}>0$, 可以看出式(31)中的可行集合是一个无界集合, 可能不存在$\pmb u_k$达到下确界. 因此, 为了方便求解, 令
$$ \begin{align} \underline{\delta}^{[0i]}\le {{\hat{\delta }}^{[0i]}} \le \bar{\delta }^{[0i]} \end{align} $$ (38) 其中, $\underline{{\delta }}^{[0i]}>0$为最小的分离阈值, $\bar{\delta }^{[0i]}$为最大的分离阈值, 虽然这样会引入一些保守性, 但是满足约束条件(36)和(38)的$\pmb u$ 仍可以视为所需要的辅助输入信号.
进一步, 优化问题(36)可以转化成
$$ \begin{align} &\min\, {\pmb u}_{k}^{\rm T} { R}{{\pmb u}_{k}}\nonumber\\[2mm] & \rm {s.t.} \quad \begin{cases} {\hat \delta ^{[0i]}} = \mathop {\min }\limits_{{\delta ^{[0i]}},{{\pmb \xi }^{[0i]}}}\, {{\delta }^{[0i]}} \\[1mm]\underline{{\delta}}^{[0i]} \le {{\hat{\delta }}^{[0i]}}\le \bar{\delta }^{[0i]} \\[1mm] {{\vartriangle }^{[0i]}}{{\pmb u}_{k}}= {G_{m}^{[0i]}}{{\pmb \xi }^{[0i]}}+{\pmb c_{m}^{[0i]}}\\[1mm] {{\left\| {{\pmb \xi }^{[0i]}} \right\|}_{\infty }} \le 1+{{\delta }^{[0i]}} \end{cases} \end{align} $$ (39) 注意到式(39)为双层优化问题, 不易得到最优解, 本文通过引入拉格朗日乘子使得式(39)中的优化问题转化为单优化问题.
定理4. 对于系统(1), 设计辅助输入信号$\pmb u_k$, 使其满足以下条件
$$ \begin{align} \begin{split} &\underset{{{\pmb u}_{k}}, {{\delta }^{[0i]}}, {{\pmb \xi }^{[0i]}}, {{\pmb \lambda }^{[0i]}}, \pmb \mu _{1}^{[0i]}, \pmb \mu _{2}^{[0i]}, \pmb b_{1}^{[0i]}, \pmb b_{2}^{[0i]}} {\mathop{\min }}\, \pmb u_{k}^{\rm T}{R}{{\pmb u}_{k}} \\ &\rm {s.t.}\\ &~~~~\left\{ \begin{aligned} &\underline{{\delta }}^{[0i]} \le {{\delta }^{[0i]}}\le \bar{\delta }^{[0i]}, \quad i=1, 2, \cdots, q \\ &{{\Delta }^{[0i]}}{{\pmb u}_{k}}={G_{m}^{[0i]}}{{\pmb \xi }^{[0i]}}+{\pmb c_{m}^{[0i]}} \\ &{{\left\| {{\pmb \xi }^{[0i]}} \right\|}_{\infty }}\le 1+{{\delta }^{[0i]}} \\ &{{(\pmb \mu _{1}^{[0i]}+\pmb \mu _{2}^{[0i]})}^{\rm T}}\pmb {1}=1 \\ &\pmb \mu _{1}^{[0i]}-\pmb \mu _{2}^{[0i]}={({G_{m}^{[0i]}})^{\rm T}}{{\pmb \lambda }^{[0i]}} \\ &b_{1, l}^{[0i]}, b_{2, l}^{[0i]}\in \{0, 1\}, \quad l=1, 2, \cdots , 2r \\ %\begin{matrix} &(\xi _{l}^{[0i]}-1-{{\delta }^{[0i]}})\in [-2(1+\bar{\delta }^{[0i]})(1-b_{1, l}^{[0i]}), 0] \\ &(\xi _{l}^{[0i]}+1+{{\delta }^{[0i]}})\in [0, 2(1+\bar{\delta }^{[0i]})(1-b_{1, l}^{[0i]})] \end{aligned} \right. \end{split} \end{align} $$ (40) 则引入辅助输入信号$u_k$ 后, 正常模型和故障模型状态所在的集合交集为空. 其中, $\pmb \lambda$, ${\pmb \mu_{1}^{[0i]}}\in{{{\bf{R}}}^{2r}}$, ${\pmb \mu_{2}^{[0i]}}$ $\in$ ${{{\bf{R}}}^{2r}}$ 为拉格朗日乘子, $b_{1, l}^{[0i]}, b_{2, l}^{[0i]}$ 为二进制变量, 向量$\pmb{1}=[\!\!\text{ }\begin{matrix} 1 & 1 & \cdots & 1 \\ \end{matrix}{{\text{ }\!\!]\!\!\text{ }}^{\rm T}}$.
证明.设拉格朗日乘子分别为${{\pmb \lambda }^{[0i]}}$, $\mu _{1, 1}^{[0i]}$, $\mu _{1, 2}^{[0i]}$, $\cdots$, $\mu _{1, 2r}^{[0i]}$, $\mu _{2, 1}^{[0i]}$, $\mu _{2, 2}^{[0i]}$, $\cdots$, $\mu _{2, 2r}^{[0i]}$.
$$ \pmb \mu _{1}^{[0i]}={\left[ \begin{matrix} \mu _{1, 1}^{[0i]} & \mu _{1, 2}^{[0i]} & \cdots & \mu _{1, 2r}^{[0i]} \end{matrix} \right]^{\rm T}} $$ (41) $$ \pmb \mu _{2}^{[0i]}={\left[ \begin{matrix} \mu _{2, 1}^{[0i]} & \mu _{2, 2}^{[0i]} & \cdots & \mu _{2, 2r}^{[0i]} \end{matrix} \right]^{\rm T}} $$ (42) 其中,
$$ \begin{align} \mu_{1, 2r}^{[0i]}\ge 0~, ~\mu_{2, 2r}^{[0i]}\ge 0 \end{align} $$ (43) 构造拉格朗日函数为
$$ \begin{align} Lag=&\ {{\delta }^{[0i]}}+{{\pmb \lambda }^{[0i]}}\left({G_{m}^{[0i]}}{{\xi }^{[0i]}}+{\pmb c_{m}^{[0i]}}-{{\vartriangle }^{[0i]}}{{\pmb u}_{k}}\right)+\notag\\ &\ \mu _{1, 1}^{[0i]}\left(-\xi _{1}^{[0i]}-1-{{\delta }^{[0i]}}\right)+\notag\\ &\ \mu _{1, 2}^{[0i]}\left(-\xi _{2}^{[0i]}-1-{{\delta }^{[0i]}}\right) +\cdots\, +\notag\\ &\ \mu _{1, 2r}^{[0i]}\left(-\xi _{2r}^{[0i]}-1-{{\delta }^{[0i]}}\right)+\notag\\ &\ \mu _{2, 1}^{[0i]}\left(\xi _{1}^{[0i]}-1-{{\delta }^{[0i]}}\right)+\notag\\ &\ \mu _{2, 2}^{[0i]}\left(\xi _{2}^{[0i]}-1-{{\delta }^{[0i]}}\right)+\notag\\ &\ \mu _{2, 2r}^{[0i]}\left(\xi _{2r}^{[0i]}-1-{{\delta }^{[0i]}}\right)\Rightarrow \notag\\ &\ {{\delta }^{[0i]}}+{{\pmb \lambda }^{[0i]}}\left({G_{m}^{[0i]}}{{\pmb \xi }^{[0i]}}+{\pmb c_{m}^{[0i]}}-\vartriangle ^{[0i]}{{\pmb u}_{k}}\right)-\notag \\ &\ {\left( \pmb \mu _{1}^{[0i]} \right)^{\rm T}}{{\pmb \xi }^{[0i]}}+ {\left(\pmb \mu _{2}^{[0i]} \right)^{\rm T}}{{\pmb \xi }^{[0i]}}- {\left( \pmb \mu _{1}^{[0i]} \right)^{\rm T}}\pmb{1}\, -\notag\\ &\ {\left( \pmb \mu _{2}^{[0i]} \right)^{\rm T}}\pmb{1}-{\left(\pmb \mu _{1}^{[0i]} \right)^{\rm T}}\pmb{1}{{\delta }^{[0i]}}-{\left( \pmb \mu _{2}^{[0i]} \right)^{\rm T}}\pmb{1}{{\delta }^{[0i]}} \end{align} $$ (44) 为了求得${{\delta }^{[0i]}}$的局部最小值, 利用文献[31]中的定理$3.3.1$和$3.4.1$, 通过对拉格朗日函数局部求导, 可以分别得到
$$ \frac{\partial Lag}{\partial {{\delta }^{[0i]}}}=1-{\left( \pmb \mu _{1}^{[0i]}+\pmb \mu _{2}^{[0i]} \right)^{\rm T}}\pmb{1}=0 $$ (45) $$ \frac{\partial Lag}{\partial {{\pmb \xi }^{[0i]}}}={\left( {G_{m}^{[0i]}} \right)^{\rm T}}{{\pmb \lambda }^{[0i]}}-\pmb \mu _{1}^{[0i]}+\pmb \mu _{2}^{[0i]}=0 $$ (46) $$ \mu _{1, l}^{[0i]}(\xi _{l}^{[0i]}-1-{{\delta }^{[0i]}})=0 , \quad l=1, 2, \cdots , 2r $$ (47) $$ \mu _{2, l}^{[0i]}(\xi _{l}^{[0i]}+1+{{\delta }^{[0i]}})=0 $$ (48) 化简式(45)和式(46), 可得:
$$ {\left( \pmb \mu _{1}^{[0i]}+\pmb \mu _{2}^{[0i]} \right)^{\rm T}}\pmb {1}=1 $$ (49) $$ \pmb \mu _{1}^{[0i]}-\pmb \mu _{2}^{[0i]}={({G_{m}^{[0i]}})^{\rm T}}{{\pmb \lambda }^{[0i]}} $$ (50) 由于式(47)和式(48)是非凸函数, 所以不能作为优化问题的约束函数. 因此, 通过引入二进制变量${\pmb b}_1^{[0i]}$, ${\pmb b}_2^{[0i]}\in {\{0, 1\}}$, 则式(47)和式(48)可以转化为
$$ \begin{align} &b_{1, l}^{[0i]}=1\Rightarrow ~~\mu _{1, l}^{[0i]}~\text{为任意值}, ~ (\xi_{l}^{[0i]}-1-{{\delta }^{[0i]}})=0 \notag \\ &b_{1, l}^{[0i]}=0\Rightarrow ~~\mu _{1, l}^{[0i]}=0, ~(\xi _{l}^{[0i]}-1-{{\delta }^{[0i]}})~\text{为任意值}\notag \\ &b_{2, l}^{[0i]}=1\Rightarrow ~~\mu _{2, l}^{[0i]}~\text{为任意值}, ~ (\xi_{l}^{[0i]}+1+{{\delta }^{[0i]}})=0 \notag \\ &b_{2, l}^{[0i]}=0\Rightarrow ~~\mu _{2, l}^{[0i]}=0, ~ (\xi_{l}^{[0i]}+1+{{\delta }^{[0i]}})~\text{为任意值} \end{align} $$ (51) 由式(43)和式(49), 可以得到$\mu _{1, l}^{[0i]}$, $\mu _{2, l}^{[0i]}\in [0, 1]$. 由式（35）可得:
$$ \begin{align} {{\left\| {{\pmb \xi }^{[0i]}} \right\|}_{\infty }}=\left\| \xi _{\max }^{[0i]} \right\| \end{align} $$ (52) 结合式(52), 式(35)可以化为
$$ \begin{align} -1-{{\bar{\delta }}^{[0i]}}\le \xi _{\max }^{[0i]}\le 1+{{\bar{\delta }}^{[0i]}} \end{align} $$ (53) 由式(47)和式(48)可得, 当$\mu _{1, l}^{[0i]}=0$ 时, $\xi _{l}^{[0i]}$ $-1-{{\delta }^{[0i]}}$ 为任意值. 由式(53) 可得, 对于任意的$\xi _{l}^{[0i]}$, 上式都成立, 因此
$$ {{\left( \xi _{l}^{[0i]} \right)}_{\max }}=1+{{\bar{\delta }}^{[0i]}} $$ (54) $$ {{\left( \xi _{l}^{[0i]} \right)}_{\min }}=-1-{{\bar{\delta }}^{[0i]}} $$ (55) 可以得到
$$ \begin{align} -2(1+\bar{\delta })\le \xi _{l}^{[0i]}-1-{{\delta }^{[0i]}}\le 0 \end{align} $$ (56) 同理可得
$$ \begin{align} 0\le \xi _{l}^{[0i]}+1+{{\delta }^{[0i]}}\le 2(1+\bar{\delta }) \end{align} $$ (57) 因此, 整理式(56)和式(57)可以得到以下边界条件
$$ \begin{align} & \mu _{1, l}^{[0i]}~, ~\mu _{2, l}^{[0i]}\in [0, 1] \nonumber\\ & (\xi _{l}^{[0i]}-1-{{\delta }^{[0i]}})\in [-2(1+{{{\bar{\delta }}}^{[0i]}}), 0] \nonumber\\ & (\xi _{l}^{[0i]}+1+{{\delta }^{[0i]}})\in [0, 2(1+{{{\bar{\delta }}}^{[0i]}})] \end{align} $$ (58) 通过式(58)中的边界条件, 以及在式(51)中引入二进制, 可以将式(47)和式(48)转化为
$$ \begin{align} \begin{split} & \mu _{1, l}^{[0i]}\le b_{1, l}^{[0i]}~, ~\mu _{2, l}^{[0i]}\le b_{2, l}^{[0i]} \\ & (\xi _{l}^{[0i]}-1-{{\delta }^{[0i]}})\in [-2(1+\bar{\delta }^{[0i]})(1-b_{1, l}^{[0i]}), 0] \\ & (\xi _{l}^{[0i]}+1+{{\delta }^{[0i]}})\in [0, 2(1+\bar{\delta }^{[0i]})(1-b_{2, l}^{[0i]})] \end{split} \end{align} $$ (59) 整理式(39), (49), (50), (59), 可得优化问题式(40).
利用正常模型和故障模型交集为空的思想设计得到辅助输入信号以后, 全对称多胞形的故障检测逻辑主要的目标是判断实际系统估计得到的状态$\hat{\pmb x}_k$ 是否发生故障, 则判断系统是否发生故障可以转化成以下逻辑问题:
$$ \begin{align} f_{{\rm result}}= \begin{cases} \text{0}, &\mbox{若}~\hat{\pmb x}_k \in {\mathcal{Z}}_k^{[0]}\\ \text{1}, &\mbox{若}~\hat{\pmb x}_k \notin {\mathcal{Z}}_k^{[0]} \end{cases} \end{align} $$ (60) 其中, $f_{{\rm result}}=0$时, 系统处于正常情况下, 否则, 系统发生故障. 式(60) 中的逻辑判断条件根据全对称多胞形的定义可以转换为以下约束问题:
$$ \begin{align} \begin{cases} {{c}_{k, 1}^{[0]}}+{{G}_{k, 1}^{[0]}}\alpha ={{\hat x}_{k, 1}} \\[2mm] {{c}_{k, 2}^{[0]}}+{{G}_{k, 2}^{[0]}}\alpha ={{\hat x}_{k, 2}} \\[1mm] ~~~~~~~~~~~~~~~~~~~~\vdots \\[1mm] {{c}_{k, n_x}^{[0]}}+{{G}_{k, n_x}^{[0]}}\alpha ={{\hat x}_{k, n_x}} \end{cases} \end{align} $$ (61) 其中, $\alpha = [\begin{matrix} {{\alpha }_{1}} & {{\alpha }_{2}} & \cdots & {{\alpha }_{r}} \\ \end{matrix}\text{ }\!]\!\!\text{ } ^{{\rm T}}$, $\| {{\alpha }_{r}} \|\le 1$. $c_{k, n_x}^{[0]}$, $G_{k, n_x}^{[0]}$分别代表全对称多胞形的中心矩阵和生成矩阵的第$n$ 行, ${\hat x}_{k, n_x}$ 代表${\hat {\pmb x}}_{k}$ 的第$n_x$ 行. 若${\hat x}_{k, n_x}$ 满足式(61)的约束条件, 则说明$\hat {\pmb x}_k$ 在全对称多胞形${\mathcal{Z}}_k^{[0]}$中, 实际系统处在正常情况下, 否则, 说明实际系统发生故障.
4. 仿真验证
为了验证本文所提方法的正确性与有效性, 采用文献[25]中的永磁直流电动机的低频线性模型进行仿真. 该系统模型为
$$ \begin{align} &\left[ \begin{matrix} \frac{{\rm d}i(t)}{{\rm d}t} \\ \frac{{\rm d}n(t)}{{\rm d}t} \\ \end{matrix} \right]=\left[ \begin{matrix} -\frac{{{R}_{a}}}{L} & -\frac{{{K}_{e}}}{L} \\ \frac{{{K}_{t}}}{{{J}_{1}}} & -\frac{{{f}_{r}}}{{{J}_{1}}} \\ \end{matrix} \right]\left[ \begin{matrix} i(t) \\ n(t) \\ \end{matrix} \right]+\left[ \begin{matrix} -\frac{1}{L} \\ 0 \\ \end{matrix} \right]u(t)\nonumber\\[2mm] &\ \pmb y(t)=\left[ \begin{matrix} 1 & 0 \\ 0 & 1 \\ \end{matrix} \right]\left[ \begin{matrix} x_1(t) \\ x_2(t) \\ \end{matrix} \right] \end{align} $$ (62) 其中, $u$, $i$, $R_a$, $L$, $K_e$, $K_t$, $J_1$, $f_r$ 分别为电枢电压, 电流, 电阻, 电感, 转矩常数, 反电动势常数, 电机惯量和摩擦系数, 转矩常数$K_t=1.0005\times K_e$. $u=$ $u_c + u_a$, 其中, $u_c$ 为系统的控制输入, $u_a$ 为系统的辅助输入. 电机处在正常情况下系统参数分别为${{R}_{a}}$ $=$ $1.2030 \Omega$, $L=5.5840\times 10^{-3}\, {\rm H}$, ${{K}_{e}}=8.1876$ $\times$ $10^{-2}\, {\rm V}$ $\cdot$ ${\rm rad}/{\rm s}$, ${{J}_{1}}=1.3528\times 10^{-4} {\rm J}$ $\cdot$ ${\rm s}^{2}/{\rm rad}$, ${{f}_{r}}=$ $2.3396\times 10^{-4} {\rm J}$ $\cdot$ ${\rm s}/{\rm rad}$. 当电机电感发生变化时, 即$L$ $=8.7548\times 10^{-3} {\rm H}$ 时, 系统发生故障. 根据电动机的机理过程, 为将电动机的转速维持在$70.3$ rad/s, 给定系统的控制输入$u_c=6$ V. 利用前向欧拉差分方法将电动机模型离散化, 其中采样时间为5 ms. 由于离散过程会使得参数引入不确定性, 因此离散化的状态动态过程增加$E{{\pmb w}_{k}}$, 其中
$$ E=\left[ \begin{matrix} -0.0254 & -0.0778 \\ -0.3996 & 0.3026 \\ \end{matrix} \right] $$ 假设系统的初始状态、测量噪声和过程干扰所在的全对称多胞形分别满足
$$ \begin{align} &\pmb x_0 \in {{\mathcal{X}}_{0}}=\left\langle \left[ \begin{matrix} 0.6 \\ 70 \\ \end{matrix} \right]\left[ \begin{matrix} 0.06 & 0 \\ 0 & 0.6 \\ \end{matrix} \right] \right\rangle\notag\\ &\pmb v_k \in{\mathcal{V}}=\left\langle \left[ \begin{matrix} 0 \\ 0 \\ \end{matrix} \right]\left[ \begin{matrix} 0.06 & 0 \\ 0 & 0.6 \\ \end{matrix} \right] \right\rangle\notag\\ &\pmb w_k \in \mathcal{W} = \langle {\bf{0}}, I_2\rangle \end{align} $$ (63) 系统的初始辅助输入信号给定为$u_a=0$. 通过求解最优化问题(40), 可得辅助输入信号, 在求解优化问题(40) 时, 取最小分离阈值${\underset{\scriptscriptstyle-}{\delta }^{[01]}=0.01}$, 最大分离阈值$\bar{\delta }^{[01]}=1.2$.
给定实际系统, 假设实际系统前50时刻发生故障的情况如下:
$$ \begin{align} \begin{cases} ~~0\le k<20, &\text{正常} \\ 20\le k<40, &\text{故障}\\ 40\le k<50, &\text{正常} \end{cases} \end{align} $$ (64) 为了验证所提方法的有效性, 本文与基于输出集合的辅助输入信号方法进行对比.
图 3是两种方法得到的辅助输入信号. 其中实线为基于本文所提方法所得的辅助输入信号, 虚线为基于输出集合交集所得到的辅助输入信号. 由图 3可以较为明显地看出, 在系统稳定后本文提出的方法得到的辅助输入信号的幅值相对较小.
图 4给出了不同方法得到的辅助输入信号对系统输出影响的对比. 其中实线表示系统未加辅助输入信号时系统的状态以及相应的输出, 虚线表示文中所提方法得到的辅助输入信号对系统状态和输出影响, 点画线表示利用输出集合得到的辅助输入信号对系统状态和输出的影响. 可以看出本文所提出的方法对原系统的影响较小.
图 4 加入辅助信号和未加辅助信号对系统输出的影响Fig. 4 Effect of adding auxiliary signal and unassisted signal on system output图 5分别为$k=20, 21, 22, 23$ 时刻以及$k=$ $40$, $41, 42, 43$ 时刻利用文中所提方法得到的辅助输入信号, 将该辅助输入信号注入系统后, 通过观测器观测得到的正常模型和故障模型的状态所在的全对称多胞形的分离情况. 其中, 三角符号以及方形符号分别为相应全对称多胞形的中心, 圆圈为实际系统观测得到的状态. 根据式(64)可得, 系统在$k= 20$ 时刻发生故障, 由图 5可以看出, 实际系统在$k=$ $21$ 时刻落在故障全对称多胞形中, 可以判断系统在$k$ $=21$ 处在故障的情况下. 同理看得, 实际系统在$k$ $=41$ 时刻落在正常全对称多胞形中. 可以判断系统在$k=41$ 处在正常的情况下. 由于实际系统在$k$=20时刻已经发生故障, 所以系统由正常到故障不会突变, 都有个缓慢的过程, 所以在$k=20$ 时刻时未检测出故障, 故障检测结果有一个时刻的延迟. 图 6为利用全对称多胞形故障检测逻辑得到的故障检测结果, 可以直观地看出利用文中方法对实际系统的故障检测结果.
图 5 状态多胞形: 通过状态交集为空的方法Fig. 5 Zonotope of the system state: the method of emptying the state intersection为了更好地展现出利用文中所提方法比基于输出集合交集为空设计辅助输入信号的保守性小, 给出了$k=20, ~21, ~22, ~23$ 相应的全对称多胞形, 见仿真图 7和图 8. 图 7和图 8为基于输出集合交集为空的主动故障检测方法所得到的结果. 图 7为系统输出所在的全对称多胞形的分离情况, 图 8为系统输出所在全对称多胞形刚好相切的情况下, 对应的状态所在全对称多胞形的情况. 可以看出正常模型和故障模型的输出所在的全对称多胞形相切时, 此时状态所在的全对称多胞形已经相离. 因此, 这说明从系统状态的角度出发, 设计辅助输入信号去激励系统实现故障分离, 该方法要比从系统输出角度出发, 设计得到的辅助信号小, 保守性小.
图 7 输出多胞形: 通过输出交集为空的方法Fig. 7 Zonotope of the system output: the method of emptying the output intersection
图 8 状态多胞形: 通过输出交集为空的方法Fig. 8 zonotope of the system state: the method of emptying the output intersection5. 结论
本文在系统过程干扰和测量噪声未知但有界的前提下, 设计了一种基于状态估计的辅助输入信号方法. 首先利用全对称多胞形卡尔曼滤波设计系统的状态观测器, 对系统的状态进行估计, 得到实时状态所在的集合. 之后在对系统影响最小的前提下, 利用正常模型的状态所在的多胞形和故障模型的状态所在的多胞形的交集为空, 设计得到实时的最优辅助输入信号. 将该辅助输入信号注入系统, 充分激励系统, 提高系统的故障检测质量. 最后通过分析数值模型验证文章所提方法. 本文所用到的设计状态观测器的方法是通过实时最小化状态所在多胞形的$F$ 半径, 该方法虽然可以实时得到状态观测器增益, 但在进行状态估计时仍不能十分准确地估计状态, 因此如何选择状态观测器较为准确估计状态仍需进一步进行研究.
-
图 1 在正常模型和故障模型下, 加入辅助信号与未加辅助信号结果的对比
Fig. 1 Comparison of the results of adding auxiliary signals and unassisted signals under normal and fault models
图 2 在正常模型和故障模型下, 加入辅助信号对状态和输出的影响
Fig. 2 The effect of adding an auxiliary signal on state and output under normal and fault models
图 4 加入辅助信号和未加辅助信号对系统输出的影响
Fig. 4 Effect of adding auxiliary signal and unassisted signal on system output
图 5 状态多胞形: 通过状态交集为空的方法
Fig. 5 Zonotope of the system state: the method of emptying the state intersection
图 7 输出多胞形: 通过输出交集为空的方法
Fig. 7 Zonotope of the system output: the method of emptying the output intersection
-
[1] Precup R E, Angelov P, Costa B S J, Sayed-Mouchaweh M. An overview on fault diagnosis and nature-inspired optimal control of industrial process applications. Computers in Industry, 2015, 74: 75-94 doi: 10.1016/j.compind.2015.03.001 [2] Monroy I, Escudero G, Graells M. Anomaly detection in batch chemical processes. Computer Aided Chemical Engineering, 2009, 26(9): 255-260 http://www.sciencedirect.com/science/article/pii/S1570794609700434 [3] Zhou M, Wang Z H, Shen Y, Shen M Q. $H_-$/${{H}}_{\infty}$ fault detection observer design in finite frequency domain for Lipschitz non-linear systems. IET Control Theory and Applications, 2017, 11(14): 2361-2369 %zhoumeng doi: 10.1049/iet-cta.2017.0357 [4] Wang Z H, Shi P, Lim C C. $H_-$/${ {H}}_{\infty}$ fault detection observer in finite frequency domain for linear parameter-varying descriptor systems. Automatica, 2017, 86: 38-45 %wangzhenhua doi: 10.1016/j.automatica.2017.08.021 [5] Han W X, Wang Z H, Shen Y, Liu Y P. ${ {H}}_{\infty}$ fault detection for linear discrete-time descriptor systems. IET Control Theory and Applications, 2018, 12(15): 2156-2163 doi: 10.1049/iet-cta.2017.1408 [6] Liu G J, Cao Y Y, Chang X H. Fault detection observer design for fuzzy systems with local nonlinear models via fuzzy Lyapunov function. International Journal of Control, Automation and Systems, 2017, 15(5): 2233-2242 doi: 10.1007/s12555-016-0496-6 [7] Zhang Q, Geng S C. Dynamic uncertain causality graph applied to dynamic fault diagnoses of large and complex systems. IEEE Transactions on Reliability, 2015, 64(3): 910- 927 doi: 10.1109/TR.2015.2416332 [8] 袁杰, 王福利, 王姝, 赵露平. 基于D-S融合的混合专家知识系统故障诊断方法. 自动化学报, 2017, 43(9): 1580-1587 doi: 10.16383/j.aas.2017.c160676Yuan Jie, Wang Fu-Li, Wang Shu, Zhao Lu-Ping. A fault diagnosis approach by D-S fusion theory and hybrid expet knowledge system. Acta Automatica Sinica, 2017, 43(9): 1580-1587 doi: 10.16383/j.aas.2017.c160676 [9] 文成林, 吕菲亚, 包哲静, 刘妹琴. 基于数据驱动的微小故障诊断方法综述. 自动化学报, 2016, 42(9): 1285-1299 doi: 10.16383/j.aas.2016.c160105Wen Cheng-Lin, Lv Fei-Ya, Bao Zhe-Jing, Liu Mei-Qin. A review of data driven-based incipient fault diagnosis. Acta Automatica Sinica, 2016, 42(9): 1285-1299 doi: 10.16383/j.aas.2016.c160105 [10] Cai B P, Zhao Y B, Liu H L, Xie M. A data-driven fault diagnosis methodology in three-phase inverters for PMSM drive systems. IEEE Transactions on Power Electronics, 2017, 32(7): 5590-5600 doi: 10.1109/TPEL.2016.2608842 [11] 刘强, 卓洁, 郎自强, 秦泗钊. 数据驱动的工业过程运行监控与自优化研究展望. 自动化学报, 2018, 44(11): 26-38 doi: 10.3969/j.issn.1003-8930.2018.11.005Liu Qiang, Zhou Jie, Lang Zi-Qiang, Qin S. Joe. Perspectives on data-driven operation monitoring and self-optimization of industrial processes. Acta Automatica Sinica, 2018, 44(11): 26-38 doi: 10.3969/j.issn.1003-8930.2018.11.005 [12] 李娟, 周东华, 司小胜, 陈茂银, 徐春红. 微小故障诊断方法综述. 控制理论与应用, 2012, 29(12): 1517-1529 https://www.cnki.com.cn/Article/CJFDTOTAL-KZLY201212000.htmLi Juan, Zhou Dong-Hua, Si Xiao-Sheng, Chen Mao-Yin, Xu Chun-Hong. Review of incipient fault diagnosis methods. Control Theory and Applications, 2012, 29(12): 1517-1529 https://www.cnki.com.cn/Article/CJFDTOTAL-KZLY201212000.htm [13] Simandal M, Puncochar I. Active fault detection and control: unified formulation and optimal design. Automatica, 2009, 45(9): 2052-2059 doi: 10.1016/j.automatica.2009.04.028 [14] Ashari A E, Nikoukhah R, Campbell S L. Auxiliary signal design for robust active fault detection of linear discrete-time systems. Automatica, 2011, 47(9): 1887-1895 doi: 10.1016/j.automatica.2011.06.009 [15] Puncochar I, Skach J. A survey of active fault diagnosis methods. IFAC Papers OnLine, 2018, 51(54): 1091-1098 [16] Wang J, Zhang J J, Qu B, Wu H Y, Zhou J L. Unified architecture of active fault detection and partial active fault-tolerant control for incipient faults. IEEE Transactions on Systems, Man, and Cybernetics Systems , 2017, 47(7): 1688 -1699 doi: 10.1109/TSMC.2017.2667683 [17] Niemann H, Poulsen N K. Active fault detection in MIMO systems. In: Proceedings of the 2014 IEEE American Control Conference. Portland, Oregon, USA: IEEE, 2014. 1975- 1980 [18] Niemann H, Poulsen N K. Active fault diagnosis in closed-loop uncertain systems. IFAC Proceedings Volumes, 2006, 39(13): 587-592 doi: 10.3182/20060829-4-CN-2909.00097 [19] Mesbah A, Streif S, Findeisen R, Braatz R D. Active fault diagnosis for nonlinear systems with probabilistic uncertainties. IFAC Proceedings Volumes, 2014, 47(3): 7079-7084 doi: 10.3182/20140824-6-ZA-1003.01594 [20] Bertsekas D P, Rhodes I B. Recursive state estimation for a set-membership description of uncertainty. IEEE Transactions on Automatic Control, 1971, 16(2): 117-128 doi: 10.1109/TAC.1971.1099674 [21] 汤文涛, 王振华, 王烨, 沈毅. 基于未知输入集员滤波器的不确定系统故障诊断. 自动化学报, 2018, 44(9): 1717-1724 doi: 10.16383/j.aas.2017.c170123Tang Wen-Tao, Wang Zhen-Hua, Wang Ye, Shen Yi. Fault diagnosis for uncertain dystems based on unknown input set-membership filters. Acta Automatica Sinica, 2018, 44(9): 1717-1724 doi: 10.16383/j.aas.2017.c170123 [22] Campbell S L, Horton K G, Nikoukhah R. Auxiliary signal design for rapid multi-model identification using optimization. Automatica, 2002, 38(8): 1313-1325 doi: 10.1016/S0005-1098(02)00040-7 [23] Zhai S C, Wang W, Ye H. Auxiliary signal design for active fault detection based on set-membership. IFAC Papers OnLine, 2015, 48(21): 452-457 doi: 10.1016/j.ifacol.2015.09.568 [24] Raimondo D M, Marseglia G R, Braatz R D, Scott J K. Closed-loop input design for guaranteed fault diagnosis using set-valued observers, Automatica, 2016, 74: 107-117 doi: 10.1016/j.automatica.2016.07.033 [25] Scott J K, Findeisen R, Braatz R D, Raimondo D M. Input design for guaranteed fault diagnosis using zonotopes, Automatica, 2014, 50: 1580-1589 doi: 10.1016/j.automatica.2014.03.016 [26] Scott J K, Findeisen R, Braatz R D, Raimondo D M. Design of active inputs for set-based fault diagnosis. In: Proceedings of the 2013 IEEE American Control Conference. Washington, DC, USA: IEEE, 2013. 3561-3566 [27] Kuhn W. Rigorously computed orbits of dynamical systems without the wrapping effect. Computing, 1998, 61(1): 47- 67 doi: 10.1007/BF02684450 [28] Wang Y, Puig V. Zonotopic extended Kalman filter and fault detection of discrete-time nonlinear systems applied to a quadrotor helicopter. In: Proceedings of the 2016 IEEE Control and Fault-Tolerant Systems, Barcelona, Spain: IEEE, 2016. 367-372 [29] Combastel C. A state bounding observer for uncertain non-linear continuous-time systems based on zonotopes. In: Proceedings of the 2005 IEEE Conference on Decision and Control and the European Control Conference, Seville, Spain: IEEE, 2005. 7228-7234 [30] Combastel C. Zonotopes and Kalman observers: Gain optimality under distinct uncertainty paradigms and robust convergence. Automatica, 2015, 55: 265-273 doi: 10.1016/j.automatica.2015.03.008 [31] Bertsekas D P. Nonlinear Programming. Belmont, MA: Athena Scientific, 1999. 314-360 期刊类型引用(1)
1. 马磊，张莹，杨春雨，王国庆. 测量值异常下多时间尺度系统H_∞滤波. 系统科学与数学. 2022(07): 1648-1659 . 
百度学术其他类型引用(14)
-
下载:
下载:
计量
- 文章访问数: 1940
- HTML全文浏览量: 162
- PDF下载量: 327
- 被引次数: 15
