Asynchronous Updating Reinforcement Learning Algorithm for Decision-making Operational Indices of Uncertain Industrial Processes
-
摘要: 运行指标决策问题是实现工业过程运行安全和生产指标优化的关键. 考虑到多运行指标决策问题求解的复杂性和工业过程生产条件动态波动引发生产指标状态的不确定性, 提出了一种策略异步更新强化学习算法自学习决策运行指标, 并给出算法收敛性的理论证明. 该算法在随机自适应动态规划框架下, 利用样本均值代替计算生产指标状态转移概率矩阵, 因此无需要求生产指标状态转移概率矩阵已知. 并且通过引入时钟和定义其阈值, 采用集中式策略评估、多策略异步更新方式用以简化求解多运行指标决策问题, 提高强化学习的学习效率. 利用可测量数据, 自学习得到的运行指标能够保证生产指标优化, 并且限制在规定范围之内. 最后, 采用中国西部某大型选矿厂的实际数据进行仿真验证, 表明该方法的有效性.Abstract: The decision-making operational index has been a key issue for achieving safe and optimal operation of industrial processes. Considering the complexity of decision making of multiple operational indices and the uncertainty of production indices caused by changes of working condition in industrial processes, this paper proposes a reinforcement learning algorithm with policy asynchronous updating for the first time aiming at self-learning operational indices, followed by the theoretical proof of convergence of the proposed algorithm. To this end, under the framework of stochastic adaptive dynamic programming, the sample mean is utilized rather than calculating the state transition probability matrix of production indices, with the outcome that the state transition probability matrix of production indices is not required to be known a priori. Distinctly from traditional synchronized policy updating, the centralized policy evaluation and asynchronous updating of multiple policies are implemented in the proposed algorithm based on the introduction of a time clock with its threshold, such that solving the concerned decision-making problem of multiple operational indices becomes easier and the learning efficiency of reinforcement learning is improved. Thus, the self-learned operational indices using measured data can ensure the optimality of production indices and limit them within the prescribed range. Experiments are conducted using the real date collected from a large-scale mineral processing plant in west China in order to illustrate the effectiveness of the approach.
-
现代控制系统变得越来越复杂, 相应地, 控制系统的安全性和可靠性就显得尤为重要. 但是, 实际系统运行时发生的各种故障会降低其可靠性甚至破坏系统的稳定性从而造成严重的安全事故. 因此, 为了提高控制系统的可靠性和安全性, 故障诊断技术在过去30年中得到了国内外学者的广泛关注, 取得了很多优秀的研究成果[1-5]. 一般来说, 故障诊断技术主要分为3个部分: 检测、隔离和估计. 故障检测的目的是确定控制系统中是否发生故障; 故障隔离则是在检测到故障后确定其位置与类型; 而故障估计则是用于得到故障的幅值信息.
故障检测作为故障诊断的第一步, 对于提高控制系统的可靠性和安全性具有十分重要的意义, 因此得到了广泛的研究. 在目前的故障检测方法中, 基于模型的故障检测方法是一种较为深入的研究方法, 涌现了很多研究成果[6-9]. 基于模型故障检测方法的主要思想是借助系统的数学模型和输入输出值设计一个故障检测观测器, 该观测器可以生成残差信号, 用于指示是否发生故障. 从其思想可以看出, 该方法的检测性能主要取决于系统数学模型和输入输出信号的准确性. 但是, 实际控制系统中总是存在着各种各样的不确定因素, 如模型不确定性、过程扰动和测量噪声等, 这些因素会降低基于模型方法的故障检测性能. 为了解决这一问题, 近年来很多学者开始研究基于模型的鲁棒故障检测方法[10-12]. 基于模型的鲁棒故障检测方法要求观测器生成的残差能够同时满足对扰动与噪声的鲁棒性和对故障的敏感性[13]. 为了实现这一性能, 基于$ H_- / H_{\infty} $的故障检测观测器设计受到了研究者们的广泛关注. 特别地, 文献[14]中给出了一种基于迭代线性矩阵不等式和多目标优化的$ H_- / H_{\infty} $故障检测观测器设计算法. 针对离散时间Takagi-Sugeno模糊系统, 文献[15]通过使用广义Kalman-Yakubovich-Popov引理设计了一类有限频$ H_- / H_{\infty} $模糊故障检测滤波器, 可以实现对传感器和执行器故障的检测. 文献[16]中设计了一种$ H_- / H_{\infty} $观测器来检测存在参数不确定性的非线性系统的故障. 针对非线性奇异系统, 文献[17]设计出了一种新型的有限频$ H_- / H_{\infty} $故障检测观测器结构, 该结构具有更多的设计自由度, 从而提高了故障检测的性能. 需要说明的是, 上述方法全都是使用$ H_{\infty} $范数来实现残差对扰动与噪声的鲁棒性. 注意到$ H_{\infty} $范数是用于描述信号能量之间的增益关系, 而实际信号的能量往往都是无界的, 故使用$ H_{\infty} $范数来描述残差的鲁棒性并不十分适合[18]. 考虑到信号虽然能量是无界的但都满足峰值有界, 因此采用描述信号峰值之间增益关系的$ L_{\infty} $范数来设计残差对扰动与噪声的鲁棒性更加合理. 特别地, 文献[19]设计了一种$ H_- / L_{\infty} $观测器来检测参数时变系统的执行器故障. 针对Lipschitz非线性系统, 文献[20]设计了一类有限频域上的$ H_- / L_{\infty} $故障检测观测器. 通过使用$ L_{\infty} $范数分析技术, 文献[21]提出了一种基于区间观测器的故障检测方法. 虽然$ L_{\infty} $范数在故障检测领域有所应用, 但与$ H_{\infty} $范数相比, 相关成果并不多. 同时, 目前大部分方法都使用$ H_- $因子来衡量残差对故障的敏感性, 同样需要考虑故障能量有界, 这一假设很难在实际应用中得到满足, 需要进行改进.
另一方面, 在故障检测观测器设计完成以后, 还需要对生成的残差进行评价以生成合理的阈值来进行故障检测. 但是, 目前的大部分方法都只选用固定的常值阈值来检测故障, 若常值阈值选取不当则容易造成虚警现象[22-24]. 虽然文献[19-20]通过使用$ L_{\infty} $范数分析技术得到了动态的故障检测阈值, 可以有效消除虚警现象, 但是其存在阈值初值过大的问题, 这会在一定程度上影响检测性能[25]. 为了解决上述问题, 学者们开始研究基于集员估计技术的动态阈值设计方法, 这类方法的主要思想是基于集员估计技术, 利用特殊的几何体如平行多面体、中心对称多面体、椭球等来包含系统每一时刻的所有无故障残差以得到动态的残差阈值集合; 然后, 通过判断实际残差是否被该阈值集合包含即可实现故障检测. 这类方法不存在虚警问题且检测性能优于文献[19-20]中的结果. 近年来, 中心对称多面体由于其求取线性映射和闵科夫斯基之和的简便性受到了学者们的关注, 涌现了一类利用中心对称多面体描述残差阈值集合进而实现故障检测的方法. 基于自适应观测器和中心对称多面体技术, 文献[26]设计了针对非线性系统的执行器故障检测方法. 文献[27]基于中心对称多面体提出了一种输入设计方法, 并将其应用到参数时变系统的故障检测中. 文献[28]在多目标观测器设计的基础上利用中心对称多面体和区间分析技术解决了线性时不变系统的传感器故障检测问题. 文献[29]将$ H_- / H_{\infty} $故障检测观测器与中心对称多面体残差评价结合来检测控制系统的执行器故障. 针对受到未知扰动和噪声影响的事件触发系统, 文献[30]设计了一种基于$ H_- / L_{\infty} $观测器与中心对称多面体残差评价的故障检测方法. 虽然上述方法利用中心对称多面体实现了较好的故障检测性能, 但中心对称多面体在其应用过程中维数会不断增长, 容易造成维数灾难. 同时, 这些方法均需要计算高维的中心对称多面体来保证故障检测精度, 计算量消耗巨大. 尽管这些计算量可以通过降维操作进行降低, 但也会在一定程度上降低检测性能, 难以很好地兼顾检测精度和计算效率[31].
除了中心对称多面体, 椭球也是集员估计方法中常用的几何体, 在控制系统状态估计领域受到了学者们的广泛关注与研究[32-37]. 特别地, 近些年来有学者借助其光滑边界的几何性质和凸优化技术研究基于椭球的故障检测方法. 针对线性变参数系统, 文献[38]通过判断包含上一时刻系统状态的椭球集与测量数据生成的超平面是否相交来检测系统故障. 通过引入适应度函数和辅助信号, 文献[39]设计了一种基于残差椭球凸优化技术的主动式故障检测方法. 针对离散时间线性切换系统, 文献[40]通过给定的输入输出数据来计算无故障情况下的系统参数椭球集, 并判断该椭球集与切换子系统参数集合是否相交来检测故障. 文献[41]通过在线求解线性矩阵不等式确定状态预测椭球和校正椭球并判断两者的交集是否为空集, 实现了对网络控制系统的故障检测. 虽然上述文献中的方法能实现较好的故障检测性能, 但大都依赖于实时求解凸优化问题或线性矩阵不等式, 这需要消耗很大的计算量, 不利于实际系统的线上实现. 同时, 判断椭球与超平面或椭球与椭球是否存在交集需要求解多组代数方程, 在一定程度上会降低故障检测方法的运算效率. 实际上, 除了利用凸优化技术或线性矩阵不等式方法来进行椭球故障检测, 还可以借助基于线性映射与闵科夫斯基之和运算的椭球分析技术来确定残差阈值集合以进行故障检测. 特别地, 椭球的线性映射与闵科夫斯基之和计算过程可以转换为简单的矩阵运算且计算过程中椭球维数不会发生改变, 无需计算高维椭球来提高精度, 计算量需求少[42]. 考虑到基于检测观测器设计和中心对称多面体残差评价的方法难以良好兼顾检测精度和计算效率, 而现有的椭球故障检测方法又大都采用凸优化或线性矩阵不等式技术, 运算效率较低, 本文旨在结合故障检测观测器设计和基于线性映射与闵科夫斯基之和运算的椭球残差评价技术, 设计一种能良好兼顾性能和计算效率的故障检测方法.
基于上述讨论, 本文针对具有未知扰动与测量噪声的线性离散时间系统, 提出了一种基于极点配置和椭球分析的传感器故障检测方法. 首先, 将传感器故障视为增广状态, 将原始系统转化为一个等效的新线性系统. 然后, 利用极点配置和$ L_{\infty} $技术设计故障检测观测器, 使得生成的残差满足故障敏感性和未知扰动鲁棒性. 最后, 使用椭球分析方法来评价残差, 生成动态的无故障残差椭球以实现故障检测. 本文的研究方法采用了基于极点配置的故障敏感性评价因子和基于$ L_{\infty} $范数的鲁棒性分析技术, 将传统$ H_- / H_{\infty} $技术需要信号能量有界的假设放松至信号峰值有界, 具有更好的适用性. 此外, 基于椭球的残差评价方法只需要简单的矩阵运算, 计算量需求小, 具备良好的计算效率, 易于实际系统实现.
1. 符号说明与相关概念
符号说明. $ {\bf{R}}^n $和$ {\bf{R}}^{n \times m} $分别表示$ n $维欧氏空间及$ n \times m $维矩阵构成的集合. $ I_n $表示$ n \times n $维的单位矩阵, $ {\boldsymbol{0}} $表示具有适当维数的零向量或零矩阵. 对于给定矩阵$ X \in {\bf{R}}^{n \times m} $, $ X^{\dagger} $代表矩阵$ X $的伪逆. 对于给定矩阵$ Y \in {\bf{R}}^{n \times n} $, $ Y^{-1} $, $ {\rm{tr}}(Y) $和$ {\rm{det}}(Y) $分别代表矩阵$ Y $的逆、迹以及行列式. 对于对称矩阵$ P\in {\bf{R}}^{n \times n} $, $ \, P \succ {\boldsymbol{0}} \; (P\prec{\boldsymbol{0}}) $表示矩阵$ P $为正定(负定)矩阵.
本文中还将使用如下的定义、性质和引理.
定义 1. 对于离散信号$ {\boldsymbol{x}} \in{\bf{R}}^{n} $, 其$ L_{\infty} $范数定义为
$$ \|{\boldsymbol{x}}\|_{\infty} = \sup\limits_{k \geq 0}\|{\boldsymbol{x}}_k\|, \ \|{\boldsymbol{x}}_k\| = {\sqrt{{\boldsymbol{x}}^{{\rm{T}}}_k{\boldsymbol{x}}_k}} $$ 定义 2. 对于两个集合$ {\boldsymbol{M}} $和$ {\boldsymbol{N}} $, 它们的闵科夫斯基和运算定义为
$$ {\boldsymbol{M}} \oplus {\boldsymbol{N}} = \left\{ m + n : m \in {\boldsymbol{M}}, n \in {\boldsymbol{N}} \right\} $$ (1) 其中, $ \oplus $表示闵科夫斯基和运算符号.
定义 3[43] . 一个$ n $维空间中的非退化椭球$ {\boldsymbol{E}}({\boldsymbol{c}} , X) $定义为
$$ {\boldsymbol{E}}({\boldsymbol{c}},X) =\{ { {\boldsymbol{x}} : ({\boldsymbol{x}}-{\boldsymbol{c}})^{{\rm{T}}}X ^{-1}({\boldsymbol{x}}-{\boldsymbol{c}}) \leq 1 }\} $$ (2) 其中, 向量$ {\boldsymbol{c}} \in{\bf{R}}^n $称为椭球$ {\boldsymbol{E}}({\boldsymbol{c}},X) $的中心, 矩阵$ X \succ {\boldsymbol{0}} \in {\bf{R}}^{n\times n} $称为椭球$ {\boldsymbol{E}}({\boldsymbol{c}},X) $的形状矩阵, 决定了椭球$ {\boldsymbol{E}}({\boldsymbol{c}},X) $的形状和体积.
注 1. 虽然式(2)中椭球的定义得到广泛使用, 但是该定义只能描述非退化形式的椭球, 即$ X \succ 0 $的情形. 事实上, 退化形式的椭球在很多实际问题中也应当被考虑到. 为了不失椭球描述的一般性, 本文引入了如下的椭球定义.
定义 4[44]. 一个$ n $维空间中的椭球$ {\cal{E}}({\boldsymbol{c}},M) $可视为一个$ n $维空间中单位球的线性映射, 定义如下:
$$ {\cal{E}}({\boldsymbol{c}},M) = \lbrace {\boldsymbol{x}} : {\boldsymbol{x}} = {\boldsymbol{c}} + M{\boldsymbol{z}}, {\boldsymbol{z}} \in{\bf{R}}^n, {\boldsymbol{z}}^{{\rm{T}}}{\boldsymbol{z}} \leq 1 \rbrace $$ (3) 其中, $ {\boldsymbol{c}} \in{\bf{R}}^n $称为椭球$ {\cal{E}}({\boldsymbol{c}},M) $的中心, $ M \in{\bf{R}}^{n\times n} $称为椭球$ {\cal{E}}({\boldsymbol{c}},M) $的形状矩阵.
注 2. 注意到式(3)中椭球的定义不要求形状矩阵$ M\succ 0 $, 故可同时描述退化与非退化形式的椭球, 具有更好的描述普适性. 特别地, 在描述非退化形式的椭球时, 定义3与定义4是等价的且$ X = MM^{{\rm{T}}} $关系成立. 不失一般性, 本文中的椭球采用式(3)中的定义形式.
性质 1[44]. 对于椭球$ {\cal{E}}({\boldsymbol{c}},M) $中的每个元素$ {\boldsymbol{x}} $进行$ {\boldsymbol{y}} = K{\boldsymbol{x}}+{\boldsymbol{b}} $的线性映射可得到一个新椭球, 且有
$$ K{\cal{E}}({\boldsymbol{c}},M) + {\boldsymbol{b}} = {\cal{E}}(K{\boldsymbol{c}}+{\boldsymbol{b}},KM)$$ (4) 其中, $ {\boldsymbol{b}}\in {\bf{R}}^n $和$ K \in {\bf{R}}^{m\times n} $为已知的向量和矩阵.
性质 2[44]. 对于给定的$ N $个椭球$ {\cal{E}}({\boldsymbol{c}}_i,M_i) $, $ i = 1, \cdots, N, $ 它们的闵科夫斯基和运算满足
$$ {\cal{E}}({\boldsymbol{c}}_1,M_1) \oplus \cdots \oplus {\cal{E}}({\boldsymbol{c}}_N,M_N) \subseteq {\cal{E}}({\boldsymbol{c}}_m,M({\boldsymbol{\alpha}})) $$ (5) 其中,
$$ {\boldsymbol{c}}_m = \sum\limits_{i = 1}^{N}{\boldsymbol{c}}_i, \ M({\boldsymbol{\alpha}})M^{{\rm{T}}}({\boldsymbol{\alpha}}) = \sum\limits_{i = 1}^{N}\frac{1}{\alpha_i}M_iM_i^{{\rm{T}}} $$ 其中, $ {\boldsymbol{\alpha}} = [\alpha_1, \, \cdots, \, \alpha_N] $ 且满足$ \sum\nolimits_{i = 1}^{N}\alpha_i = 1 $, $ \alpha_i >0 $, $ \forall \, i = 1, \cdots, N $.
注 3. 不同的参数$ \alpha_1, \, \cdots, \, \alpha_N $会产生不同形状和体积的椭球$ {\cal{E}}({\boldsymbol{c}}_m, M({\boldsymbol{\alpha}})) $. 为了保证椭球闵科夫斯基之和运算的精度, 可选择不同的优化准则来求取参数$ \alpha_1, \cdots, \alpha_N ,$ 以得到优化的椭球$ {\cal{E}}({\boldsymbol{c}}_m,M({\boldsymbol{\alpha}})) $. 常用的优化准则为椭球的体积和半轴长平方和, 分别对应最小化椭球形状矩阵的行列式和迹. 在这两种优化准则下, 优化参数$ \alpha_1, \, \cdots, \, \alpha_N $具体形式为[44]:
1) 体积最小化椭球$ {\cal{E}}({\boldsymbol{c}}_m,M({\boldsymbol{\alpha}}^+)) $的优化参数向量$ {\boldsymbol{\alpha}}^+ = [\alpha_1^+, \cdots, \alpha_N^+] $可由下式得到
$$ \alpha_i^+ = \arg\min\limits_{0 <\alpha_i^+ < 1}{\rm{lndet}}(\hat{M}_i\hat{M}_i^{{\rm{T}}}), \,\; \hat{M}_1 = M_1 $$ (6) $$ \hat{M}_{i+1}\hat{M}_{i+1}^{{\rm{T}}} = {\frac{1}{\alpha_i^+}}\hat{M}_i\hat{M}_i^{{\rm{T}}}+ {\frac{1}{1-\alpha_i^+}}M_{i+1}M_{i+1}^{{\rm{T}}} $$ (7) 2) 半轴长平方和最小化椭球$ {\cal{E}}({\boldsymbol{c}}_m,M({\boldsymbol{\alpha}}^*)) $的优化参数向量$ {\boldsymbol{\alpha}}^* = [\alpha_1^*, \cdots, \alpha_N^*] $可由下式得到:
$$ \alpha_i^* = \frac{\sqrt{{\rm{tr}}(M_iM_i^{{\rm{T}}})}}{\sum\limits_{i = 1}^{N}\sqrt{{\rm{tr}}(M_iM_i^{{\rm{T}}})}}, \ \; i = 1, \cdots, N $$ (8) 从式(6) ~ (8)中可以看出, 求解最小化体积椭球的参数过程十分复杂, 耗费很大计算量, 不利于实际系统实现. 与之相比, 求取最小化半轴长平方之和椭球的参数计算量较小, 实时性较好, 因此本文选用椭球的半轴长平方之和作为优化准则来求取最优的闵科夫斯基之和椭球$ {\cal{E}}({\boldsymbol{c}}_m,M({\boldsymbol{\alpha}})) $.
引理 1[45]. 对于矩阵$ X\in{\bf{R}}^{a\times b}, \ Y\in $ ${\bf{R}}^{b\times c} , \ Z\in{\bf{R}}^{a\times c} $. 如果$ {\rm{rank}}(Y) = c $, 则方程$ XY = Z $的通解为
$$ X = ZY^{\dagger} + S[I_b-YY^{\dagger}] $$ (9) 其中, $ S\in{\bf{R}}^{a\times b} $为任意矩阵.
2. 问题描述
考虑如下存在传感器故障的线性离散时间系统
$$ \left\{\begin{aligned} &{\boldsymbol{x}}_{k+1} = A{\boldsymbol{x}}_k + B{\boldsymbol{u}}_k+D_w{\boldsymbol{w}}_k\\ &{\boldsymbol{y}}_k = C{\boldsymbol{x}}_k+F{\boldsymbol{f}}_k+D_v{\boldsymbol{v}}_k \end{aligned}\right. $$ (10) 其中, $ {\boldsymbol{x}}_k\in{\bf{R}}^{n_x} $, $ {\boldsymbol{u}}_k \in {\bf{R}}^{n_u} $, $ {\boldsymbol{y}}_k \in {\bf{R}} ^{n_y} $分别是系统的状态向量、控制输入和测量输出, $ {\boldsymbol{f}}_k\in{\bf{R}}^{n_f} $表示传感器故障, $ {\boldsymbol{w}}_k\in{\bf{R}}^{n_w} $和$ {\boldsymbol{v}}_k \in{\bf{R}}^{n_v} $分别为系统受到的未知过程干扰和测量噪声; $ A, B, C, D_w, D_v, F $是具有适当维数的常数矩阵. 不失一般性, 假设矩阵$ F $为列满秩的, 即$ {\rm{rank}}(F) = n_f $.
本文假设系统(10)的状态变量初值和受到的过程干扰与测量噪声均为未知但有界的, 即
$$ \left\{\begin{aligned} &\|{\boldsymbol{x}}_0-{\boldsymbol{c}}_0\| \leq \tilde{x}_0 \\ &\|{\boldsymbol{w}}_k\| \leq \tilde{w} = \|{\boldsymbol{w}}\|_{\infty} \\ &\|{\boldsymbol{v}}_k\| \leq \tilde{v} = \|{\boldsymbol{v}}\|_{\infty} \end{aligned}\right. $$ (11) 其中, $ {\boldsymbol{c}}_0 \in {\bf{R}}^{n_x} $为已知向量, $\tilde{x} _0$, $ \tilde{w} $和$ \tilde{v} $为已知常数.
为了检测故障$ {\boldsymbol{f}}_k $, 本文将其视为增广状态, 则可得到如下的增广状态向量
$$ \bar{{\boldsymbol{x}}}_k = \begin{bmatrix} {\boldsymbol{x}}_k \\ {\boldsymbol{f}}_k \end{bmatrix} $$ (12) 并构造出如下的增广系统
$$ \left\{\begin{aligned} &\bar{{\boldsymbol{x}}}_{k+1} = \bar{A}\bar{{\boldsymbol{x}}}_k + \bar{B}{\boldsymbol{u}}_k+\bar{D}_w{\boldsymbol{w}}_k+\bar{F}{\boldsymbol{f}}_{k+1}\\ &{\boldsymbol{y}}_k = \bar{C}\bar{{\boldsymbol{x}}}_k+D_v{\boldsymbol{v}}_k \end{aligned}\right. $$ (13) 其中,
$$ \begin{array}{l} \bar{A} = \begin{bmatrix} A &{\boldsymbol{ 0}} \\ {\boldsymbol{0}} & {\boldsymbol{0 }} \end{bmatrix},\, \bar{B} = \begin{bmatrix} B \\ {\boldsymbol{0}} \end{bmatrix}\\ \bar{C} = \begin{bmatrix} C & F \end{bmatrix},\, \bar{D}_w = \begin{bmatrix} D_w \\ {\boldsymbol{0}} \end{bmatrix},\, \bar{F} = \begin{bmatrix} {\boldsymbol{0}} \\ I_{n_f} \end{bmatrix}\end{array} $$ 显然, 上述状态增广过程并未采用任何假设, 所以系统(13)与原系统(10)完全等价. 因此, 若系统(13)检测到故障, 则表明原系统(10)发生了故障. 所以, 原系统(10)的传感器故障检测问题就转变为了对增广系统(13)的故障检测.
为了实现上述目标, 本文针对系统 (13)提出了一种结合极点配置和椭球分析技术的故障检测方法. 首先, 针对系统 (13)设计了一个故障检测观测器, 其产生的残差能够同时满足对故障的敏感性和对扰动与噪声的鲁棒性. 特别地, 残差对故障的敏感性通过极点配置方法来实现; 同时使用$ L_{\infty} $技术来抑制扰动与噪声对残差的影响. 基于设计的检测观测器, 本文采用椭球技术来进行残差评价并给出故障检测结果.
3. 主要成果
3.1 故障检测观测器设计
针对系统(13)构造如下形式的故障检测观测器
$$ \left\{ \begin{array}{l} \hat{\bar{{\boldsymbol{x}}}}_{k+1} = \bar{A}\hat{\bar{{\boldsymbol{x}}}}_k+\bar{B}{\boldsymbol{u}}_k+L({\boldsymbol{y}}_k-\bar{C}\hat{\bar{{\boldsymbol{x}}}}_k) \\ {\boldsymbol{r}}_k = {\boldsymbol{y}}_k -\bar{C}\hat{\bar{{\boldsymbol{x}}}}_k \end{array} \right. $$ (14) 其中, $ \hat{\bar{{\boldsymbol{x}}}}_k \in {\bf{R}}^{n_x+n_f} $是增广状态的估计向量, $ {\boldsymbol{r}}_k \in {\bf{R}}^{n_y} $为故障检测残差, $ L\in{\bf{R}}^{(n_x+n_f)\times n_y} $是待设计的故障检测观测器增益矩阵.
定义估计误差为
$$ {\boldsymbol{e}}_k = \bar{{\boldsymbol{x}}}_k - \hat{\bar{{\boldsymbol{x}}}}_k $$ (15) 则由式(13)和式(14), 可得如下的误差系统
$$ \left\{\begin{aligned} &{\boldsymbol{e}}_{k+1} = (\bar{A}-L\bar{C}){\boldsymbol{e}}_k+\bar{F}{\boldsymbol{f}}_{k+1}\;+ \\ &\qquad\quad\;\bar{D}_w{\boldsymbol{w}}_k-LD_v{\boldsymbol{v}}_k \\ &{\boldsymbol{r}}_k = \bar{C}{\boldsymbol{e}}_k+D_v{\boldsymbol{v}}_k \end{aligned}\right. $$ (16) 为了分析传感器故障和扰动与噪声对残差的影响, 将误差动态系统 (16) 拆分为如下的两个子系统
$$ \left\{\begin{aligned} &{\boldsymbol{e}}_{f,k+1} = (\bar{A}-L\bar{C}){\boldsymbol{e}}_{f,k}+\bar{F}{\boldsymbol{f}}_{k+1}\\ &{\boldsymbol{r}}_{f,k} = \bar{C}{\boldsymbol{e}}_{f,k} \end{aligned}\right. $$ (17) $$ \left\{\begin{aligned} &{\boldsymbol{e}}_{d,k+1} = (\bar{A}-L\bar{C}){\boldsymbol{e}}_{d,k}+\bar{D}_w{\boldsymbol{w}}_k\; -\\ & \qquad \qquad LD_v{\boldsymbol{v}}_k\\ &{\boldsymbol{r}}_{d,k} = \bar{C}{\boldsymbol{e}}_{d,k}+D_v{\boldsymbol{v}}_k \end{aligned}\right. $$ (18) 其中, $ {\boldsymbol{e}}_{d,0} = {\boldsymbol{e}}_0 $, $ {\boldsymbol{e}}_{f,0} = 0 $, 则有$ {\boldsymbol{e}}_k = {\boldsymbol{e}}_{d,k}+{\boldsymbol{e}}_{f,k} $和$ {\boldsymbol{r}}_k = {\boldsymbol{r}}_{d,k}+{\boldsymbol{r}}_{f,k} $.
假设存在一个常数$\zeta$使得下式成立
$$ (\bar{A}-L\bar{C})\bar{F} = \zeta \bar{F} $$ (19) 从式(19)可看出, $ \zeta $为矩阵$ \bar{A}-L\bar{C} $的特征值, 也为误差系统(16)的极点. 为了保证误差系统(16)的稳定性, 常数$ \zeta $需在区间$ (0,1) $中选择, 即 $ 0< \zeta<1 $.
将式(19)代入到误差系统(17)中并进行迭代, 可得到
$$ {\boldsymbol{e}}_{f,k} = \zeta^{k-1}\bar{C}\bar{F}{\boldsymbol{f}}_1 + \cdots + \zeta \bar{C}\bar{F}{\boldsymbol{f}}_{k-1} + \bar{C}\bar{F}{\boldsymbol{f}}_{k}$$ (20) 注意到$ F = \bar{C}\bar{F} $, 则式(20) 等价为
$$ {\boldsymbol{e}}_{f,k} = \zeta^{k-1}F{\boldsymbol{f}}_1+\cdots + \zeta F{\boldsymbol{f}}_{k-1} +F{\boldsymbol{f}}_{k} $$ (21) 从式(21)中可看出, $ \zeta $反映了残差对故障的敏感程度, 可通过提高$ \zeta $的数值来提升残差对故障的敏感性.
为了求取矩阵$ L $, 将式(19)改写为
$$ \bar{A}\bar{F}-\zeta \bar{F} = L\bar{C}\bar{F} $$ (22) 注意到$ \bar{C}\bar{F} = F $且矩阵$ F $满足列满秩条件, 则根据引理1可知满足式(19)的矩阵$ L $的通解为
$$ L = \Theta_1 + S\Theta_2 $$ (23) 其中, $ S \in {\bf{R}}^{(n_x+n_f) \times n_y} $为任意选择的常数矩阵, $ \Theta_1 $和$ \Theta_2 $的具体形式为
$$\left\{\begin{aligned} &\Theta_1 = (\bar{A}\bar{F}-\zeta \bar{F})(\bar{C}\bar{F})^{\dagger} \\ &\Theta_2 = I_{n_y}-(\bar{C}\bar{F})(\bar{C}\bar{F})^{\dagger} \end{aligned}\right.$$ 注意到式(23)仅给出了增益矩阵$ L $满足故障敏感性能的通解形式, 并没有给出矩阵$ S $的选取方法. 同时, 故障检测过程中还需要使用鲁棒技术来抑制未知扰动与噪声对残差的影响. 因此, 本文将引入$ L_{\infty} $技术来保证误差系统(18)是渐近稳定的且残差$ {\boldsymbol{r}}_{d,k} $对扰动与噪声满足如下的$ L_{\infty} $性能
$$ \|{\boldsymbol{r}}_{d,k}\| \leq \sqrt{(\gamma_w+\gamma_v)(\lambda(1-\lambda)^k V_0 + \gamma_w \tilde{w}^2+ \gamma_v \tilde{v}^2}) $$ (24) 其中, $ \gamma_w>0, $ $ \gamma_v>0 ,$ $ 0 < \lambda < 1 $为给定的常数, $ V_0 = {\boldsymbol{e}}^{{\rm{T}}}_0P{\boldsymbol{e}}_0 $, $ P \succ {\boldsymbol{0}} \in {\bf{R}}^{(n_x+n_f)\times (n_x+n_f)} $为待设计的正定矩阵.
基于误差系统(17)和(18), 本文给出定理1, 保证残差满足式(21)和式(24)中的故障敏感性能和未知扰动与噪声鲁棒性能.
定理 1. 对于给定常数$ \gamma_w>0, \gamma_v>0 ,0 <\lambda < $$ 1 $和$ 0<\zeta<1, $ 如果存在常数$ \mu>0 ,$ 正定矩阵$ P \in {\bf{R}}^{(n_x+n_f)\times (n_x+n_f)} $和矩阵$ W \in {\bf{R}}^{(n_x+n_f)\times n_f}, $ 使得如下不等式成立
$$ \left[\begin{array}{*{20}{c}} (\lambda-1)P & 0 & 0 & \Omega_1 \\ \star & -\mu I_{n_w} & 0 & \Omega_2 \\ \star & \star & -\mu I_{n_v} & \Omega_3 \\ \star & \star & \star & -P \end{array}\right] \prec {\boldsymbol{0}} $$ (25) $$ \left[\begin{array}{*{20}{c}} \lambda P & 0 & 0 & \bar{C}^{{\rm{T}}} \\ \star & \Gamma_w & 0 & 0\\ \star & \star & \Gamma_v & D_v^{{\rm{T}}} \\ \star & \star & \star & \Gamma_y \end{array}\right] \succ {\boldsymbol{0}} $$ (26) 其中,
$$ \begin{array}{l} \Omega_1 = (P\bar{A}-P\Theta_1\bar{C}-W\Theta_2\bar{C})^{{\rm{T}}}, \ \Omega_2 = (P\bar{D}_w)^{{\rm{T}}} \\ \Omega_3 = -(P\Theta_1D_v+W\Theta_2D_v)^{{\rm{T}}},\ \Gamma_w = (\gamma_w - \mu)I_{n_w} \\ \Gamma_v = (\gamma_v - \mu)I_{n_v}, \ \Gamma_y = (\gamma_w+\gamma_v) I_{n_y} \end{array} $$ 则残差$ {\boldsymbol{r}}_k $满足式(21)和式(24)中的故障敏感性能和未知扰动鲁棒性能. 特别地, 当线性矩阵不等式(25)和(26)可解时, 则矩阵$ S $和$ L $可由下式求得
$$ S = P^{-1}W, \ L = \Theta_1 + S\Theta_2 $$ (27) 证明. 选取如下形式的李雅普洛夫函数
$$ V_k = {\boldsymbol{e}}_{d,k}^{{\rm{T}}}P{\boldsymbol{e}}_{d,k}, \ P \succ {\boldsymbol{0}} $$ (28) 则根据误差系统(18), 可得
$$ \begin{split} & \Delta V_k = V_k -V_{k+1} = {\boldsymbol{e}}^{{\rm{T}}}_{d,k+1}P{\boldsymbol{e}}_{d,k+1} -{\boldsymbol{e}}_{d,k}^{{\rm{T}}}P{\boldsymbol{e}}_{d,k} =\\ &\quad \qquad \begin{bmatrix} {\boldsymbol{e}}_{d,k} \\ {\boldsymbol{w}}_k \\ {\boldsymbol{v}}_k \end{bmatrix}^{{\rm{T}}} \begin{bmatrix} \Upsilon_{11} & \Upsilon_{12} & \Upsilon_{13} \\ \star & \Upsilon_{22} & \Upsilon_{23} \\ \star & \star & \Upsilon_{33} \end{bmatrix} \begin{bmatrix} {\boldsymbol{e}}_{d,k} \\ {\boldsymbol{w}}_k \\ {\boldsymbol{v}}_k \end{bmatrix}\\[-20pt] \end{split} $$ (29) 其中,
$$ \begin{array}{l} \Upsilon_{11} = (\bar{A}-L\bar{C})^{{\rm{T}}}P(\bar{A}-L\bar{C})- P \\ \Upsilon_{12} = (\bar{A}-L\bar{C})^{{\rm{T}}}P\bar{D}_w\\ \Upsilon_{13} = (\bar{A}-L\bar{C})^{{\rm{T}}} P(-LD_v) \\ \Upsilon_{22} = \bar{D}_w^{{\rm{T}}}P\bar{D}_w \\ \Upsilon_{23} = \bar{D}_w^{{\rm{T}}}P(-LD_v) \\ \Upsilon_{33} = (-LD_v)^{{\rm{T}}}P(-LD_v) \end{array} $$ 注意到$ W = PS $, $ L = \Theta_1 + S\Theta_2 $, 则式(25)可写为
$$ \left[ \begin{array}{*{20}{c}} (\lambda-1)P & 0 & 0 & (P(\bar{A}-L\bar{C}))^{{\rm{T}}} \\ \star & -\mu I_{n_w} & 0 & (P\bar{D}_w)^{{\rm{T}}}\\ \star & \star & -\mu I_{n_v} & (-PLD_v)^{{\rm{T}}} \\ \star &\star &\star & -P \end{array}\right] \prec {\boldsymbol{0}} $$ 在上式两边左乘右乘矩阵$ \Xi $和$ \Xi^{{\rm{T}}} $可得:
$$ \begin{bmatrix} \Upsilon_{11} & \Upsilon_{12} & \Upsilon_{13} \\ \star & \Upsilon_{22} & \Upsilon_{23} \\ \star &\star & \Upsilon_{33} \end{bmatrix} + \begin{bmatrix} \lambda P & 0 & 0 \\ \star & -\mu I_{n_w} & 0 \\ \star & \star & -\mu I_{n_v} \end{bmatrix} \prec {\boldsymbol{0}}$$ (30) 其中, 矩阵$ \Xi $具体形式为
$$ \Xi = \begin{bmatrix} I_{n_x+n_f} & 0 & 0 & (\bar{A}-L\bar{C})^{{\rm{T}}}\\ 0 & I_{n_w} & 0 & \bar{D}_w^{{\rm{T}}} \\ 0 & 0 & I_{n_v} & (-LD_v)^{{\rm{T}}} \end{bmatrix} $$ 不等式(30)两边左乘右乘向量 $ [{\boldsymbol{e}}_{d,k}^{{\rm{T}}} \quad {\boldsymbol{w}}^{{\rm{T}}}_k \quad {\boldsymbol{v}}^{{\rm{T}}}_k] $及其转置可得
$$ \Delta V_k < -\lambda V_k+\mu {\boldsymbol{w}}^{{\rm{T}}}_k{\boldsymbol{w}}_k + \mu {\boldsymbol{v}}^{{\rm{T}}}_k{\boldsymbol{v}}_k $$ (31) 当扰动$ {\boldsymbol{w}}_k $和噪声$ {\boldsymbol{v}}_k $均为零时, 由不等式(31)可推导出
$$ \Delta V_k = V_{k+1}-V_k < -\lambda V_k < 0 $$ (32) 这表明误差系统(18)是渐近稳定的.
另一方面, 不等式(31)等价于
$$ V_{k+1} < (1-\lambda)V_k+\mu {\boldsymbol{w}}^{{\rm{T}}}_k{\boldsymbol{w}}_k + \mu {\boldsymbol{v}}^{{\rm{T}}}_k{\boldsymbol{v}}_k$$ 则可得
$$ V_k < (1-\lambda)V_{k-1}+\mu \tilde{w}^2 +\mu \tilde{v}^2 $$ 通过迭代上式可得
$$ \begin{split} V_k <\ &(1-\lambda)^kV_0+\mu \sum\limits^{k-1}_{i = 0}(1-\lambda)^i(\tilde{w}^2+\tilde{v}^2) \leq\\ & (1-\lambda)^kV_0+\mu \frac{(1-\lambda)^k}{\lambda}\tilde{w}^2 + \mu \frac{(1-\lambda)^k}{\lambda}\tilde{v}^2 \leq \\ & (1-\lambda)^kV_0 + \frac{\mu \tilde{w}^2}{\lambda}+ \frac{\mu \tilde{v}^2}{\lambda}\\[-15pt] \end{split} $$ (33) 对不等式(26)运用舒尔补引理[46], 可将其转换为
$$ \begin{split} & \begin{bmatrix} \lambda P & 0 & 0 \\ \star & (\gamma_{w}-\mu)I_{n_w} & 0 \\ \star & \star & (\gamma_{v}-\mu)I_{n_v} \end{bmatrix} - \\ & \qquad\qquad \frac{1}{\gamma_{w}+\gamma_{v}} \begin{bmatrix} \bar{C}^{{\rm{T}}} \\ 0 \\ D_v^{{\rm{T}}} \end{bmatrix} \begin{bmatrix} \bar{C} & 0 & D_v \end{bmatrix} \succ 0\end{split} $$ (34) 不等式(34)两边左乘右乘向量 $ [{\boldsymbol{e}}_{d,k}^{{\rm{T}}} \quad {\boldsymbol{w}}^{{\rm{T}}}_k \quad {\boldsymbol{v}}^{{\rm{T}}}_k] $及其转置可得
$$ \begin{array}{l} \begin{bmatrix} {\boldsymbol{e}}_{d,k} \\ {\boldsymbol{w}}_k \\ {\boldsymbol{v}}_k \end{bmatrix}^{{\rm{T}}} \left( \begin{bmatrix} \lambda P & 0 & 0 \\ \star & (\gamma_{w}-\mu)I_{n_w} & 0 \\ \star & \star & (\gamma_{v}-\mu)I_{n_v} \end{bmatrix} \right. \nonumber -\\ \qquad\quad \left. \dfrac{1}{\gamma_{w}+\gamma_{v}} \begin{bmatrix} \bar{C}^{{\rm{T}}} \\ 0 \\ D_v^{{\rm{T}}} \end{bmatrix} \begin{bmatrix} \bar{C} & 0 & D_v \end{bmatrix} \right) \begin{bmatrix} {\boldsymbol{e}}_{d,k} \\ {\boldsymbol{w}}_k \\ {\boldsymbol{v}}_k \end{bmatrix} > 0 \end{array} $$ 根据式(18)和式(28), 上式等价于如下形式
$$ \begin{split} \|{\boldsymbol{r}}_{d,k}\|^2 <&\ (\gamma_{w}+\gamma_{v})(\lambda V_k + (\gamma_{w}-\mu){\boldsymbol{w}}^{{\rm{T}}}_k{\boldsymbol{w}}_k\ +\\ &(\gamma_{v}-\mu){\boldsymbol{v}}^{{\rm{T}}}_k{\boldsymbol{v}}_k) \leq\\ &(\gamma_{w}+\gamma_{v})(\lambda V_k +(\gamma_{w}-\mu)\tilde{w}^2 \ +\\ & (\gamma_{v} - \mu)\tilde{v}^2) \end{split} $$ 将式(33)代入上式可得
$$ \begin{split} \|{\boldsymbol{r}}_{d,k}\|^2 \leq&\ (\gamma_{w}+\gamma_{v})\Biggr(\lambda\left((1-\lambda)^kV_0 + \frac{\mu \tilde{w}^2}{\lambda} \right.+\\ &\left. \frac{\mu \tilde{v}^2}{\lambda}\right) \ + (\gamma_{w}-\mu)\tilde{w}^2 + (\gamma_{v} - \mu)\tilde{v}^2 \Biggr)=\\ & (\gamma_{w}+\gamma_{v})(\lambda(1-\lambda)^k V_0 + \gamma_{w}\tilde{w}^2+\gamma_{v}\tilde{v}^2) \end{split} $$ 上式等价于
$$ \|{\boldsymbol{r}}_{d,k}\| \leq \sqrt{(\gamma_w+\gamma_v)(\lambda(1-\lambda)^k V_0 + \gamma_w \tilde{w}^2+ \gamma_v \tilde{v}^2}) $$ 由此可知, 残差$ {\boldsymbol{r}}_k $满足式(24)中的$ L_{\infty} $鲁棒性能.
□ 注4. 为了使设计的故障检测观测器(14)性能最佳, 增益矩阵的设计是一个多目标优化问题. 一方面, 需要使得残差对故障的敏感性好, 即$ \zeta $的数值尽可能的大; 另一方面, 需要抑制未知扰动和噪声对残差的影响, 即$ \gamma_w $和$ \gamma_v $的数值应尽可能的小. 注意到$ 0<\zeta<1 $, 因此可在该区间内选取线性矩阵不等式(25)和(26)有解的$ \zeta $的最大值. 同时, 对于给定的$ \zeta $, 可通过求解下列优化问题来确定故障检测观测器的增益矩阵
$$ \min \ (\gamma_w+\gamma_v), \qquad {\rm{s.t.}}\ \ (25),(26) $$ (35) 若优化问题(35)可解, 则观测器(14)的增益矩阵$ L $可由$ S = P^{-1}W $, $ L = \Theta_1 + S\Theta_2 $得到.
3.2 椭球残差评价
基于模型的故障检测方法主要包括两个步骤: 残差生成和残差评价, 这两部分对于故障检测性能都有重要的影响. 但是, 现有文献中的大部分工作都集中在残差生成器的设计上, 只有少部分的工作涉及到残差评价部分. 在本节中, 本文提出了一种基于椭球分析的动态检测阈值生成方法, 可以高效快速地实现对残差的评价. 需要说明的是, 本节假设故障检测观测器已由第3.1节中提出的方法设计完成, 即增益矩阵$ L $已知并在此基础上进行椭球残差评价.
从式(17)和式(18)中可看出, 在理想状态下, 无故障时残差等于零. 然而实际系统中总是存在着未知扰动和噪声, 因此无故障时残差往往不为零, 但是可以使用椭球来包住所有无故障时残差的可能值. 首先, 根据椭球定义将有界假设(11)转化为如下的椭球形式
$$ {\boldsymbol{x}}_0 \in {\cal{E}}({\boldsymbol{c}}_0,\tilde{M}_0) , \ {\boldsymbol{w}}_k \in {\cal{E}}(0,W), \ {\boldsymbol{v}}_k \in {\cal{E}}(0,V)$$ 其中, $ \tilde{M}_0 = \tilde{x}_0I_{n_x} $, $W = \|{\boldsymbol{w}}\|_\infty I_{n_w}$, $V = \|{\boldsymbol{v}}\|_\infty I_{n_v}$.
初始时刻, 系统通常不发生故障, 即有${\boldsymbol{f}}_0 = 0,$ 则根据增广向量$ \bar{{\boldsymbol{x}}}_k $的定义可知存在如下从属关系
$$ \bar{{\boldsymbol{x}}}_0 \in {\cal{E}}(\bar{{\boldsymbol{c}}}_0,M_0) $$ 其中,
$$ \bar{{\boldsymbol{c}}}_0 = \begin{bmatrix} {\boldsymbol{c}}_0 \\ 0 \end{bmatrix}, \ M_0 = \begin{bmatrix} \tilde{M}_0 & 0 \\ 0 & 0 \end{bmatrix} $$ $ \quad \ $基于上述的椭球有界假设, 本文给出如下定理来求解包住所有无故障时残差可能值的动态椭球.
定理2. 给定初值$ \bar{{\boldsymbol{x}}}_0 \in {\cal{E}}(\bar{{\boldsymbol{c}}}_0,M_0) $和$ \hat{\bar{{\boldsymbol{x}}}}_0 = \bar{{\boldsymbol{c}}}_0, $ 则误差系统(18)任意$ k $时刻的无故障残差$ {\boldsymbol{r}}_{d,k} $可被椭球$ {\cal{E}}(0,R_k) $包含, 即有$ {\boldsymbol{r}}_{d,k} \in {\cal{E}}(0,R_k) $且形状矩阵$ R_k $满足如下迭代式
$$ M_{k+1}M_{k+1}^{{\rm{T}}} = \frac{1}{\alpha_{1,k}^*}\Pi_{x,k} + \frac{1}{\alpha_{2,k}^*}\Pi_w + \frac{1}{\alpha_{3,k}^*}\Pi_l $$ (36) $$ R_kR_k^{{\rm{T}}} = \frac{1}{\beta_{1,k}^*}\Pi_{c,k} + \frac{1}{\beta_{2,k}^*}\Pi_v \qquad\qquad \quad \qquad$$ (37) 其中,
$$ \begin{split} &\beta_{1,k}^* = \frac{\sqrt{{\rm{tr}}(\Pi_{c,k})}}{\Sigma_{1,k}}, \ \beta_{2,k}^* = \frac{\sqrt{{\rm{tr}}(\Pi_v)}}{\Sigma_{1,k}} \\ &\alpha_{1,k}^* = \frac{\sqrt{{\rm{tr}}(\Pi_{x,k})}}{\Sigma_{2,k}}, \ \alpha_{2,k}^* = \frac{\sqrt{{\rm{tr}}(\Pi_w)}}{\Sigma_{2,k}} \\ &\alpha_{3,k}^* = \frac{\sqrt{{\rm{tr}}(\Pi_l)}}{\Sigma_{2,k}}, \ \Sigma_{1,k} = \sqrt{{\rm{tr}}(\Pi_{c,k})} + \sqrt{{\rm{tr}}(\Pi_v)} \\ &\Sigma_{2,k} = \sqrt{{\rm{tr}}(\Pi_{x,k})} + \sqrt{{\rm{tr}}(\Pi_w)} + \sqrt{{\rm{tr}}(\Pi_l)} \\ &\Pi_{x,k} = (\bar{A}-L\bar{C})M_kM_k^{{\rm{T}}}(\bar{A}-L\bar{C})^{{\rm{T}}} \\ &\Pi_w = \bar{D}_w^{{\rm{T}}}WW^{{\rm{T}}}\bar{D}_w^{{\rm{T}}}, \ \Pi_l = (LD_v)VV^{{\rm{T}}}(LD_v)^{{\rm{T}}} \\ &\Pi_{c,k} = \bar{C}M_kM_k^{{\rm{T}}}\bar{C}^{{\rm{T}}}, \ \Pi_v = D_vVV^{{\rm{T}}}D_v^{{\rm{T}}} \end{split} $$ 证明. 给定初值$ \bar{{\boldsymbol{x}}}_0 \in {\cal{E}}(\bar{{\boldsymbol{c}}}_0,M_0) $和$ \hat{\bar{{\boldsymbol{x}}}}_0 = \bar{{\boldsymbol{c}}}_0 $, 则由式(15)可得
$$ {\boldsymbol{e}}_0 = \bar{{\boldsymbol{x}}}_0 - \hat{\bar{{\boldsymbol{x}}}}_0 \in {\cal{E}}(\bar{{\boldsymbol{c}}}_0,M_0) \oplus (-\hat{\bar{{\boldsymbol{x}}}}_0) = {\cal{E}}(0,M_0) $$ 因为, ${\boldsymbol{w}}_k \in {\cal{E}}(0,W),\;\;{\boldsymbol{v}}_k \in {\cal{E}}(0,V) ,\;\; {\boldsymbol{e}}_0 \in {\cal{E}}(0, M_0)$且无故障时有${\boldsymbol{e}}_{k} = {\boldsymbol{e}}_{d,k}$, 故可知存在$ {\boldsymbol{e}}_{d,k}\in{\cal{E}}(0, M_k) $. 另一方面, 根据椭球性质1可得到如下等式
$$ \begin{split} &(\bar{A}-L\bar{C}){\boldsymbol{e}}_{d,k} \in (\bar{A}-L\bar{C}) {\cal{E}}(0,M_k) =\\ & \qquad{\cal{E}}(0,(\bar{A}-L\bar{C})M_k) \\ & \bar{D}_w{\boldsymbol{w}}_k \in \bar{D}_w{\cal{E}}(0,W) = {\cal{E}}(0,\bar{D}_wW) \\ &-LD_v{\boldsymbol{v}}_k \in -LD_v{\cal{E}}(0,V) = {\cal{E}}(0,-LD_vV) \\ &\bar{C}{\boldsymbol{e}}_{d,k} \in \bar{C} {\cal{E}}(0,M_k) = {\cal{E}}(0,\bar{C}M_k) \\ &D_v{\boldsymbol{v}}_{k} \in D_v{\cal{E}}(0,V) = {\cal{E}}(0,D_vV) \end{split} $$ 另外, 通过误差动态系统(18), 可得
$$ \begin{split} &{\boldsymbol{e}}_{d,k+1} = (\bar{A}-L\bar{C}){\boldsymbol{e}}_{d,k}+\bar{D}_w{\boldsymbol{w}}_k-LD_v{\boldsymbol{v}}_k \in\\ &\qquad\qquad{\cal{E}}(0,(\bar{A}-L\bar{C})M_k) \oplus {\cal{E}}(0,\bar{D}_wW)\;\oplus\\ &\qquad \qquad {\cal{E}}(0,-LD_vV) \\ &{\boldsymbol{r}}_{d,k} = \bar{C}{\boldsymbol{e}}_{d,k}+D_v{\boldsymbol{v}}_k \in {\cal{E}}(0,\bar{C}M_k) \oplus {\cal{E}}(0,D_vV) \end{split} $$ 通过对上式使用椭球性质2, 可得
$$ \begin{split} &{\boldsymbol{e}}_{d,k+1} \in {\cal{E}}(0,(\bar{A}-L\bar{C})M_k) \oplus {\cal{E}}(0,\bar{D}_wW) \;\oplus\\ &\qquad{\cal{E}}(0,-LD_vV) \subseteq {\cal{E}}(0,M({\boldsymbol{\alpha}}_k)) \\ &{\boldsymbol{r}}_{d,k} \in {\cal{E}}(0,\bar{C}M_k) \oplus {\cal{E}}(0,D_vV) \subseteq {\cal{E}}(0,R({\boldsymbol{\beta}}_k)) \end{split} $$ 其中, $ {\boldsymbol{\alpha}}_k\, =\, [\alpha_{1,k} \quad \alpha_{2,k} \quad \alpha_{3,k}],{\boldsymbol{\beta}}_k \,=\,[\beta_{1,k} \quad \beta_{2,k}], $矩阵$ M({\boldsymbol{\alpha}}_k) $和$ R({\boldsymbol{\beta}}_k) $满足如下形式:
$$ \begin{split} &M({\boldsymbol{\alpha}}_k){M^{{\rm{T}}}({\boldsymbol{\alpha}}_k)} = \frac{1}{\alpha_{1,k}}\Pi_{x,k} + \frac{1}{\alpha_{2,k}}\Pi_w + \frac{1}{\alpha_{3,k}}\Pi_l\\ &R({\boldsymbol{\beta}}_k){R^{{\rm{T}}}({\boldsymbol{\beta}}_k)} = \frac{1}{\beta_{1,k}}\Pi_{c,k} + \frac{1}{\beta_{2,k}}\Pi_v\end{split} $$ 根据注3中的优化方法, 可得参数的最优解$ {\boldsymbol{\alpha}}_k^* $和$ {\boldsymbol{\beta}}_k^* $为
$$ \begin{split} &\alpha_{1,k}^* = \frac{\sqrt{{\rm{tr}}(\Pi_{x,k})}}{\Sigma_{2,k}}, \,\; \alpha_{2,k}^* = \frac{\sqrt{{\rm{tr}}(\Pi_w)}}{\Sigma_{2,k}} \\ &\alpha_{3,k}^* = \frac{\sqrt{{\rm{tr}}(\Pi_l)}}{\Sigma_{2,k}},\; \beta_{1,k}^* = \frac{\sqrt{{\rm{tr}}(\Pi_{c,k})}}{\Sigma_{1,k}}\;\\ &\beta_{2,k}^* = \frac{\sqrt{{\rm{tr}}(\Pi_v)}}{\Sigma_{1,k}} \end{split} $$ 最后, 令$ M_{k+1} = M({\boldsymbol{\alpha}}_k^*) $ 和 $ R_k = R({\boldsymbol{\beta}}_k^*) $, 即可得证.
□ 在系统无故障时, 存在$ {\boldsymbol{r}}_{d,k} = {\boldsymbol{r}}_{k} $关系成立, 即无故障时始终有从属关系$ {\boldsymbol{r}}_{k} \in {\cal{E}}(0,R_k) $成立. 但是, 当系统发生故障后, $ {\boldsymbol{r}}_{k} \in {\cal{E}}(0,R_k) $的条件就无法保证. 由此, 本文给出了如下的故障检测逻辑:
$$ \left\{\begin{aligned} &{\boldsymbol{r}}_{k} \in {\cal{E}}(0,R_k), \qquad \qquad \sigma_k = 0 \\ &{\boldsymbol{r}}_{k} \not \in {\cal{E}}(0,R_k), \qquad \qquad \sigma_k = 1 \end{aligned}\right. $$ (38) 其中, $ \sigma_k $为故障检测结果指示值, $ \sigma_k = 0, $系统正常; $ \sigma_k = 1, $系统故障. 特别地, 故障检测逻辑(38)可通过求解如下带约束线性规划问题来进行判断
$$ {\boldsymbol{r}}_{k} = R_k{\boldsymbol{z}}_r, \ {\boldsymbol{z}}_r \in {\bf{R}}^{n_y},\ {\boldsymbol{z}}_r^{{\rm{T}}}{\boldsymbol{z}}_r \leq 1 $$ (39) 若规划问题(39)有解, 则$ {\boldsymbol{r}}_{k} \in {\cal{E}}(0,R_k) ,$ $ \sigma_k = 0, $ 系统正常; 反之, 则$ {\boldsymbol{r}}_{k} \not \in {\cal{E}}(0,R_k), $ $ \sigma_k = 1 ,$ 系统故障.
注5. 从定理2中可看出, 动态椭球$ {\cal{E}}(0,R_k) $的计算只需要简单的矩阵运算且无需降维操作, 更不涉及到凸优化问题或线性矩阵不等式的实时求解, 具备良好的计算效率. 同时, 对无故障残差动态椭球生成方法的时间计算复杂度的比较与分析结果具体如下:
矩阵$ L ,$ $ \bar{A}-L\bar{C}, $ $ \Pi_{w}, $ $ \Pi_{l} $和$ \Pi_{v} $均为常值, 可通过线下计算提前确定, 不占用线上运行时的时间计算复杂度. 在线上运行时, 矩阵$ \Pi_{x,k} $和$ \Pi_{c,k} $的运算比其他参数的求取复杂的多, 故无故障残差值动态椭球方法的时间计算复杂度主要来源于对矩阵$ \Pi_{x,k} $和$ \Pi_{c,k} $的求取. 注意到矩阵$ M_kM_k^{\rm{T}} $和$ R_kR_k^{\rm{T}} $可当作整体来进行计算, 则求取矩阵$ \Pi_{x,k} $和$ \Pi_{c,k} $的时间计算复杂度为O($ (n_x+n_f)^3+(n_x+n_f)n_y^2 $), 故本文中残差椭球生成算法的时间计算复杂度量级为$ {\rm{O}}(n^3) $, 对计算量需求较小. 与之相比, 通过实时线上求解凸优化问题或线性矩阵不等式来确定椭球集的方法则需要消耗巨大的计算量. 以文献[41]中的方法为例, 该方法需要在线求解两个线性矩阵不等式以得到状态的预测椭球和校正椭球. 根据文献[46]中的评估准则可知, 求解这两个椭球所需的时间计算复杂度约为$ {\rm{O}}((\dfrac{(n_x+1)n_x^3}{2})^{2.5}) $和$ {\rm{O}}((\dfrac{(n_x+1)n_x^2n_y^2}{2}(1+n_v+n_x))^{2.5}) $, 计算复杂度量级很高, 计算量需求远高于本文方法. 另一方面, 本文的故障检测逻辑无需判断椭球与超平面或椭球与椭球的交集存在性, 一定程度上也降低了计算量. 因此, 基于椭球分析的残差评价方法具备良好的计算效率, 易于实际系统实现.
注6. 为了理论简便性起见, 本文只考虑了扰动和噪声所在椭球的中心为原点且形状矩阵为常值的情况, 即$ {\boldsymbol{w}}_k \in {\cal{E}}(0,W) $和$ {\boldsymbol{v}}_k \in {\cal{E}}(0,V) $. 事实上, 本文中的残差评价方法也可以推广到扰动和噪声所在椭球的中心为非原点且形状矩阵时变的情形, 即$ {\boldsymbol{w}}_k \in {\cal{E}}({\boldsymbol{c}}_{w,k},W_k) $和$ {\boldsymbol{v}}_k \in {\cal{E}}({\boldsymbol{c}}_{v,k},V_k). $ 此时, 定理2推广如下.
定理3. (定理2的推广)给定初值$ \bar{{\boldsymbol{x}}}_0 \in {\cal{E}}(\bar{{\boldsymbol{c}}}_0,M_0) $和$ \hat{\bar{{\boldsymbol{x}}}}_0 = \bar{{\boldsymbol{c}}}_0 ,$ 则误差系统(18)任意$ k $时刻的无故障残差$ {\boldsymbol{r}}_{d,k} $可被椭球$ {\cal{E}}({\boldsymbol{c}}_{r,k},R_k) $包含, 即有${\boldsymbol{r}}_{d,k} \in $${\cal{E}}({\boldsymbol{c}}_{r,k}, R_k) $且椭球中心$ {\boldsymbol{c}}_{r,k} $和形状矩阵$ R_k $满足如下的迭代等式:
$$ \begin{split} &{\boldsymbol{c}}_{e,k+1} = (\bar{A}-L\bar{C}){\boldsymbol{c}}_{e,k}+\bar{D}_w{\boldsymbol{c}}_{w,k}-LD_v{\boldsymbol{c}}_{v,k}\\ &{\boldsymbol{c}}_{r,k} = \bar{C}{\boldsymbol{c}}_{e,k}+D_v{\boldsymbol{c}}_{v,k}, \ {\boldsymbol{c}}_{e,k} = 0 \\ &M_{k+1}M_{k+1}^{{\rm{T}}} = \frac{1}{\alpha_{1,k}^*}\Pi_{x,k} + \frac{1}{\alpha_{2,k}^*}\Pi_{w,k} + \frac{1}{\alpha_{3,k}^*}\Pi_{l,k}\\ &R_kR_k^{{\rm{T}}} = \frac{1}{\beta_{1,k}^*}\Pi_{c,k} + \frac{1}{\beta_{2,k}^*}\Pi_{v,k} \end{split} $$ 其中,
$$ \begin{split} &\beta_{1,k}^* = \frac{\sqrt{{\rm{tr}}(\Pi_{c,k})}}{\Sigma_{1,k}}, \ \beta_{2,k}^* = \frac{\sqrt{{\rm{tr}}(\Pi_{v,k})}}{\Sigma_{1,k}} \\ &\alpha_{1,k}^* = \frac{\sqrt{{\rm{tr}}(\Pi_{x,k})}}{\Sigma_{2,k}}, \ \alpha_{2,k}^* = \frac{\sqrt{{\rm{tr}}(\Pi_{w,k})}}{\Sigma_{2,k}} \\ &\alpha_{3,k}^* = \frac{\sqrt{{\rm{tr}}(\Pi_{l,k})}}{\Sigma_{2,k}},\ \Sigma_{1,k} = \sqrt{{\rm{tr}}(\Pi_{c,k})} + \sqrt{{\rm{tr}}(\Pi_{v,k})} \\ &\Sigma_{2,k} = \sqrt{{\rm{tr}}(\Pi_{x,k})} + \sqrt{{\rm{tr}}(\Pi_{w,k})} + \sqrt{{\rm{tr}}(\Pi_{l,k})} \\ &\Pi_{x,k} = (\bar{A}-L\bar{C})M_kM_k^{{\rm{T}}}(\bar{A}-L\bar{C})^{{\rm{T}}} \end{split} $$ $$ \begin{split}& \Pi_{w,k} = \bar{D}_w^{{\rm{T}}}W_kW_k^{{\rm{T}}}\bar{D}_w^{{\rm{T}}}, \ \Pi_{l,k} = (LD_v)V_kV_k^{{\rm{T}}}(LD_v)^{{\rm{T}}} \\ &\Pi_{c,k} = \bar{C}M_kM_k^{{\rm{T}}}\bar{C}^{{\rm{T}}}, \ \Pi_{v,k} = D_vV_kV_k^{{\rm{T}}}D_v^{{\rm{T}}}\end{split} $$ 需要说明的是, 在未知扰动和噪声所在椭球的中心为非原点且形状矩阵时变时, 无故障残差动态椭球生成算法的计算量会略有上升, 但是时间计算复杂度量级仍为$ {\rm{O}}(n^3) $, 即所设计方法依旧保持着较好的计算效率, 利于线上实时应用.
4. 仿真结果
本节通过一个二阶RC电路模型来验证所提出传感器故障检测方法的有效性与优越性, 该二阶RC电路的原理图如图1 所示.
二阶RC电路的动态模型可描述为
$$ \begin{split} \begin{bmatrix} \dot{U}_1(t) \\ \dot{U}_2(t) \end{bmatrix} =& \begin{bmatrix} -\dfrac{R_1+R_2}{R_1R_2C_1} & \dfrac{1}{R_2C_1} \\ \dfrac{1}{R_2C_2} & -\dfrac{1}{R_2C_2} \end{bmatrix} \begin{bmatrix} U_1(t) \\ U_2(t) \end{bmatrix} +\\ &\begin{bmatrix} \dfrac{1}{R_1C_1} \\ 0 \end{bmatrix}U(t) \end{split}$$ 其中, $ U(t) $表示电路的输入电压, $ U_1(t) $和$ U_2(t) $分别表示电容$ C_1 $和$ C_2 $两端的电压.
动态模型(4)可写为如下的线性连续时间系统
$$ \left\{\begin{aligned} &\dot{{\boldsymbol{x}}}(t) = A_c{\boldsymbol{x}}(t) + B_c{\boldsymbol{u}}(t) \\ &{\boldsymbol{y}}(t) = C_c{\boldsymbol{x}}(t) \end{aligned}\right. $$ (40) 其中,
$$ \begin{split} &{\boldsymbol{x}}(t) = \begin{bmatrix} U_1(t) \\ U_2(t) \end{bmatrix}, \ {\boldsymbol{y}}(t) = \begin{bmatrix} U_1(t) \\ U_1(t)+U_2(t) \end{bmatrix}\\ &{\boldsymbol{u}}(t) = U(t), \ A_c = \begin{bmatrix} -\dfrac{R_1+R_2}{R_1R_2C_1} & \dfrac{1}{R_2C_1} \\ \dfrac{1}{R_2C_2} & -\dfrac{1}{R_2C_2} \end{bmatrix} \\ &B_c = \begin{bmatrix} \dfrac{1}{R_1C_1} \\ 0 \end{bmatrix}, \ C_c = \begin{bmatrix} 1 & 0 \\ 1 & 1 \end{bmatrix} \end{split} $$ 仿真中, 二阶RC电路的电子元件参数为$ R_1 = R_2 = 200\; \Omega $, $C_1 = C_2 = 1\;000\; \text{μ}{\rm{F}}$, 同时考虑系统受到的未知扰动和噪声以及传感器故障, 再利用欧拉一步化离散方法(离散时间$ T_s = 0.05\;{\rm{s}} $)将二阶RC电路系统(40)转换为式(10)形式的离散系统, 相关参数矩阵为
$$ \begin{array}{l} A = \begin{bmatrix} 0.50 & 0.25 \\ 0.25 & 0.75 \end{bmatrix},\ B = \begin{bmatrix} 0.25 \\ 0 \end{bmatrix},\ C = \begin{bmatrix} 1 & 0 \\ 1 & 1 \end{bmatrix} \\ F = \begin{bmatrix} 1 & 0 \\ 1 & 1 \end{bmatrix},\ D_w = 0.1I_2,\ D_v = 0.02I_2 \end{array} $$ 则可得广义系统(13)的矩阵参数为
$$ \begin{array}{l} \bar{A} = \begin{bmatrix} 0.50 & 0.25 & 0 & 0 \\ 0.25 & 0.75 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix},\; \bar{B} = \begin{bmatrix} 0.25 \\ 0 \\ 0 \\ 0 \end{bmatrix}\\ \bar{C} = \begin{bmatrix} 1 & 0 & 1 & 0 \\ 1 & 1 & 1 & 1 \end{bmatrix} , \bar{D}_w = \begin{bmatrix} 0.1 & 0 \\ 0 & 0.1 \\ 0 & 0 \\ 0 & 0 \end{bmatrix} , \bar{F} = \begin{bmatrix} 0 & 0 \\ 0 & 0 \\ 1 & 0 \\ 0 & 1 \end{bmatrix} \end{array} $$ 选取$ \zeta = 0.75 , \lambda = 0.5 $, 并求解优化问题(35), 则可得到$ \mu = 0.3021 , \gamma_w = 0.3022 ,\gamma_{v} = 0.3346 $以及观测器增益矩阵$ L $为
$$ L = \begin{bmatrix} 0.1443 & 0.1640 \\ 0.0584 & 0.2409 \\ -0.5781 & -0.0147 \\ 0.7745 & -0.4167 \end{bmatrix} $$ 仿真中, 选取系统(13)的状态初值为$\bar{{\boldsymbol{x}}}_0 = [0.1$ $ 0\;\;0 \;\; 0]^{{\rm{T}}} ,$ 控制输入为$ {\boldsymbol{u}}_k = 3\sin(0.5k), $ 未知干扰和测量噪声为$ {\boldsymbol{w}}_k = [0.2\sin(k) \quad 0.2\cos(0.5k)]^{{\rm{T}}} $和$ {\boldsymbol{v}}_k = [0.2\sin(k) \quad 0.2\cos(0.5k)]^{{\rm{T}}} .$ 此外, 选取状态估计初值为$ \hat{\bar{{\boldsymbol{x}}}}_0 = \bar{{\boldsymbol{c}}}_0 = [0 \quad 0 \quad 0 \quad 0]^{{\rm{T}}} ,$ 相关的椭球形状矩阵为$M_0 = {\rm{diag}}\{0.1,0.1,0,0\} ,$ $ W = 0.2I_2 ,$ $ V = 0.2I_2 $.
本文主要考虑两种常见类型的传感器故障: 突变故障和时变故障. 首先, 考虑如下的传感器突变故障
$$ {\boldsymbol{f}}_k = \left\{\begin{aligned} & [0 \quad 0]^{{\rm{T}}}, \qquad \quad 0 \leq k < 100\\ & [0.1 \quad 0]^{{\rm{T}}}, \qquad\, 100 \leq k \leq 200 \end{aligned}\right. $$ 传感器突变故障检测的仿真结果如图2和图3所示. 图2中给出了故障检测结果指示值$ \sigma_k $的曲线, 从图2中可看出, 在传感器发生突变故障后, 本文所设计的方法可以快速地检测到故障, 具备高效的故障检测性能. 为了直观理解基于椭球分析的残差评价方法, 图3 中给出了采样时刻$ k $在95到105的时间内, 残差$ {\boldsymbol{r}}_k $与残差椭球$ {\cal{E}}(0,R_k) $的仿真结果. 从图3中可看出, $ {\boldsymbol{r}}_{k} \not \in {\cal{E}}(0,R_k), k \geq 100 $. 根据检测逻辑(38)可知检测到系统发生故障, 验证了椭球残差评价方法的有效性.
图 3 突变故障的残差${\boldsymbol{r}}_k$与残差椭球${\cal{E}}(0,R_k)$ ($k = 95\sim105$)Fig. 3 Residual ${\boldsymbol{r}}_k$ and residual ellipsoid ${\cal{E}}(0,R_k)$ of abrupt fault ($k = 95\sim105$)仿真中, 考虑如下形式的时变传感器故障
$$ {\boldsymbol{f}}_k = \left\{\begin{aligned} & [0 \quad 0]^{{\rm{T}}}, \qquad \qquad \;\;\qquad \qquad 0 \leq k < 100 \\ & [0 \quad 0.2+0.1\sin(0.5k)]^{{\rm{T}}}, \quad \, 100 \leq k \leq 200 \end{aligned}\right. $$ 时变传感器故障检测的仿真结果如图4和图5所示. 图4中给出了故障检测结果指示值$ \sigma_k $的曲线, 从图4中可以看出, 本文所设计的方法在系统传感器发生时变故障后可以快速给出准确的检测结果, 具备良好的检测性能. 图5 中给出了采样时刻$ k $在95到105时间内, 残差$ {\boldsymbol{r}}_k $与残差椭球$ {\cal{E}}(0,R_k) $的仿真结果. 从图5中可以看出, $ {\boldsymbol{r}}_{k} \not \in {\cal{E}}(0,R_k), k \geq 100 $, 同样表明系统发生了故障, 再次验证了椭球残差评价方法的有效性.
图 5 时变故障的残差${\boldsymbol{r}}_k$与残差椭球${\cal{E}}(0,R_k)$ ($k = 95\sim105$)Fig. 5 Residual ${\boldsymbol{r}}_k$ and residual ellipsoid ${\cal{E}}(0,R_k)$ of time-varying fault ($k = 95\sim105$)为了验证本文检测方法的优越性, 仿真中将本文方法与$ H_- / H_{\infty} $和$ H_- / L_{\infty} $故障检测方法进行对比.
针对系统(10), 设计如下形式的$ H_- / H_{\infty} $故障检测观测器
$$ \left\{ \begin{array}{l} \hat{{\boldsymbol{x}}}_{k+1} = A\hat{\boldsymbol{x}}_k+B{\boldsymbol{u}}_k+L_h({\boldsymbol{y}}_k-C\hat{ {\boldsymbol{x}}}_k) \\ {\boldsymbol{\varrho}}_k = {\boldsymbol{y}}_k -C\hat{{\boldsymbol{x}}}_k \end{array} \right. $$ (41) 其中, $ \hat{{\boldsymbol{x}}}_k \in {\bf{R}}^{n_x} $是系统状态估计向量, $ {\boldsymbol{\varrho}}_k \in {\bf{R}}^{n_y} $为故障检测残差, $ L_h\in{\bf{R}}^{n_x\times n_y} $是待设计的增益矩阵. 通过使用$ H_- / H_{\infty} $设计方法可得到矩阵$ L_h $为
$$ L_h = \begin{bmatrix} 0.5639 & 0.0636 \\ -0.9843 & 1.0206 \end{bmatrix} $$ 针对系统(10), 设计如下形式的$ H_- / L_{\infty} $故障检测观测器
$$ \left\{ \begin{array}{l} \widehat{{\boldsymbol{x}}}_{k+1} = A\widehat{\boldsymbol{x}}_k+B{\boldsymbol{u}}_k+L_l({\boldsymbol{y}}_k-C\widehat{ {\boldsymbol{x}}}_k) \\ {\boldsymbol{\varsigma}}_k = {\boldsymbol{y}}_k -C\widehat{{\boldsymbol{x}}}_k \end{array} \right. $$ (42) 其中, $ \widehat{{\boldsymbol{x}}}_k \in {\bf{R}}^{n_x} $是系统状态估计向量, $ {\boldsymbol{\varsigma}}_k \in {\bf{R}}^{n_y} $为故障检测残差, $ L_l\in{\bf{R}}^{n_x\times n_y} $是待设计的增益矩阵. 通过使用$ H_- / L_{\infty} $设计方法可得到矩阵$ L_l $为
$$ L_l = \begin{bmatrix} 0.5539 & 0.0763 \\ -0.9429 & 0.9944 \end{bmatrix} $$ 仿真中, $ H_- / H_{\infty} $和$ H_- / L_{\infty} $故障检测方法也采用椭球技术来进行残差评价, 并采用如下的检测逻辑
$$ \left\{\begin{aligned} &{\boldsymbol{\varrho}}_{k} \in {\cal{E}}(0,H_k), \qquad \qquad \delta_k = 0 \\ &{\boldsymbol{\varrho}}_{k} \not \in {\cal{E}}(0,H_k), \qquad \qquad \delta_k = 1 \end{aligned}\right. $$ (43) $$ \left\{\begin{aligned} &{\boldsymbol{\varsigma}}_{k} \in {\cal{E}}(0,L_k), \qquad \qquad \, \epsilon_k = 0 \\ &{\boldsymbol{\varsigma}}_{k} \not \in {\cal{E}}(0,L_k), \qquad \qquad \, \epsilon_k = 1\end{aligned}\right. $$ (44) 其中, $ {\cal{E}}(0,H_k) $为$ H_- / H_{\infty} $故障检测方法的残差椭球, $ {\cal{E}} (0,L_k) $为$ H_- / L_{\infty} $故障检测方法的残差椭球, $ \delta_k $和$ \epsilon_k $为两种方法故障检测结果的指示值. 故障检测逻辑(43)和(44)同样通过求解带约束的线性规划问题来进行判断, 此处不再赘述.
不失公平性, 本文方法与 $ H_- / H_{\infty} $和$ H_-/ L_{\infty} $故障检测方法均采用与上述相同的仿真参数. 仿真中, 考虑如下的传感器微小突变故障:
$$ {\boldsymbol{f}}_k = \left\{\begin{aligned} & [0 \quad 0]^{{\rm{T}}}, \qquad \qquad 0 \leq k < 100 \\ & [0.03 \quad 0]^{{\rm{T}}} ,\qquad \;\; 100 \leq k \leq 200 \end{aligned}\right. $$ 仿真对比结果如图6 ~ 9所示. 图6中给出了3种方法故障检测结果指示值$ \sigma_k $、$ \delta_k $和$ \epsilon_k $的对比曲线. 从图中可看出, 本文中的方法对传感器的微小突变故障仍有检测能力, $ H_- / L_{\infty} $方法虽然能检测到故障, 但漏检率明显高于本文方法; 而$ H_- / H_{\infty} $方法则完全无法检测到该微小突变故障. 仿真结果说明本文方法具备更良好的故障检测性能. 图7 ~ 9中分别给出了采样时刻$ k $在95到105时间内, 本文方法的残差$ {\boldsymbol{r}}_k $与残差椭球$ {\cal{E}}(0,R_k) $、$ H_- / H_{\infty} $方法的残差$ {\boldsymbol{\varrho}}_k $与残差椭球$ {\cal{E}}(0,H_k) $以及$ H_- / L_{\infty} $方法的残差$ {\boldsymbol{\varsigma}}_k $与残差椭球$ {\cal{E}}(0,L_k) $的仿真结果. 从图中可看出, ${\boldsymbol{r}}_{k} \not \in {\cal{E}}(0,R_k) ,\;103 \leq k \leq 105$, 而$ {\boldsymbol{\varrho}}_{k} \in{\cal{E}}(0,H_k), $${\boldsymbol{\varsigma}}_{k} \in {\cal{E}}(0,L_k), 95 \leq k \leq 105 $. 这表明本文方法检测到了故障, 但$ H_- / H_{\infty} $和$ H_- / L_{\infty} $方法在此时间段内未检测到故障, 验证了本文方法的有效性与优越性.
图 6 微小突变故障检测指示值对比结果Fig. 6 Comparison result of small abrupt fault detection indication values
图 7 本文方法的残差${\boldsymbol{r}}_k$与残差椭球${\cal{E}}(0,R_k)$ ($k = 95\sim105$)Fig. 7 Residual ${\boldsymbol{r}}_k$ and residual ellipsoid ${\cal{E}}(0,R_k)$ by the proposed approach ($k = 95\sim105$)
图 8 $H_-/H_{\infty}$方法的残差${\boldsymbol{\varrho}}_k$与残差椭球${\cal{E}}(0,H_k)$ ($k = 95\sim105$)Fig. 8 Residual ${\boldsymbol{\varrho}}_k$ and residual ellipsoid ${\cal{E}}(0,H_k)$ by the $H_-/H_{\infty}$ method ($k = 95\sim105$)
图 9 $H_-/L_{\infty}$方法的残差${\boldsymbol{\varsigma}}_k$与残差椭球${\cal{E}}(0,L_k)$ ($k = 95\sim105$)Fig. 9 Residual ${\boldsymbol{\varsigma}}_k$ and residual ellipsoid ${\cal{E}}(0,L_k)$ by the $H_-/L_{\infty}$ method ($k = 95\sim105$)为了进一步说明本文方法在故障检测性能上的优势, 仿真中将本文方法与文献[29]中的中心多面体故障检测方法进行对比. 文献[29]中的方法设计了一个多目标故障检测观测器来生成检测残差$ {\boldsymbol{\phi}}_k $; 之后, 利用中心多面体和区间分析技术来得到无故障时残差所处的区间$ [\underline{{\boldsymbol{\phi}}}_k,\overline{{\boldsymbol{\phi}}}_k] $; 最后, 通过判断残差$ {\boldsymbol{\phi}}_k $是否在区间$ [\underline{{\boldsymbol{\phi}}}_k,\overline{{\boldsymbol{\phi}}}_k] $中来检测故障, 故障检测逻辑为
$$ \left\{\begin{aligned} &{\boldsymbol{\phi}}_k \in [\underline{{\boldsymbol{\phi}}}_k,\overline{{\boldsymbol{\phi}}}_k], \qquad \qquad \varepsilon_k = 0 \\ &{\boldsymbol{\phi}}_k \not\in [\underline{{\boldsymbol{\phi}}}_k,\overline{{\boldsymbol{\phi}}}_k], \qquad \qquad \varepsilon_k = 1\end{aligned}\right. $$ (45) 其中, $ \varepsilon_k $为文献[29]中方法故障检测结果指示值, $ \varepsilon_k = 0 $, 系统正常; $ \varepsilon_k = 1 $, 系统故障.
仿真对比过程中, 采用与上述相同的仿真参数和传感器微小突变故障, 由此可得到如图10 所示的残差区间分析仿真结果和如图11 所示的故障检测对比结果. 通过对故障检测数据统计后可知, 在传感器故障的100个采样时刻内, 本文方法共在43个采样时刻检测到故障; 文献[29]中方法共在30个采样时刻检测到故障, 说明本文中的椭球故障检测方法可以得到比文献[29]中方法更良好的故障检测结果, 漏检率更低, 验证了本文方法优秀的故障检测性能.
为了说明本文方法计算量需求小、运算效率高的优势, 仿真中将本文方法与文献[41]中的椭球故障检测方法进行对比. 对比过程中, 采用与上述相同的仿真参数和传感器微小突变故障, 将两种方法在同一环境下各运行7次并记录算法的运行时间, 由此得到如表1 所示的数据对比结果. 仿真对比实验选用了如下测试环境: Intel(R) Core(TM) i7-7700HQ CPU@ 2.8 GHz、16.00 GB 内存、Windows10系统. 从表1 中可以看出, 本文方法运行时间更短, 可得到比文献[41]中方法更高的计算效率. 这一对比结果也验证了注5中算法时间计算复杂度的分析结果, 文献[41]中高量级时间计算复杂度方法所需的运行时间远大于本文中$ {\rm{O}}(n^3) $时间计算复杂度的方法.
仿真中还对比了本文方法与文献[41]中方法的故障检测性能, 仿真对比结果如图12 所示. 通过对故障检测数据统计后可知, 在传感器微小突变故障作用的100个采样时刻内, 本文方法共在43个采样时刻检测到故障, 文献[41]中方法共在27个采样时刻检测到故障, 说明本文方法可以得到比文献[41]中方法更优秀的故障检测性能, 验证了本文方法的有效性.
5. 结束语
本文针对具有未知扰动与测量噪声的线性离散时间系统, 提出了一种新的传感器故障检测方法. 首先, 通过将传感器故障视为增广状态, 将原始系统转换为一个等效的新线性系统. 基于该系统, 本文使用鲁棒观测器设计和极点配置方法构造了一个故障检测观测器, 使得生成的残差能够同时满足对扰动与噪声的鲁棒性和对故障的敏感性, 并将检测观测器的设计问题转化为易于求解的线性矩阵不等式形式. 为了检测故障, 本文还提出了一种基于椭球分析的残差评价方法, 可通过判断残差是否被无故障残差椭球包含来检测故障. 特别地, 无故障残差椭球生成算法所需的计算量较小, 具备良好的计算效率, 易于实际系统实现. 最后, 通过一个二阶RC电路的仿真算例验证了所提出方法的有效性与优越性. 此外, 本文中所提出线性系统传感器故障检测方法可以通过线性变参数模型近似的方式推广到非线性系统, 这将是我们下一步的研究工作.
-
图 1 工业过程运行指标决策问题
Fig. 1 Decision-making problem of operational indices in industrial processes
图 3 多执行－评判结构下运行指标自学习决策流程图
Fig. 3 Flowchart of self-learning decision making of operational indices with multiple actors-critic structure
图 5 精矿产量和精矿品位损失函数
Fig. 5 Loss functions of the concentrate yield and concentrate grade
图 11 策略异步更新和策略同步更新强化学习算法时间消耗对比
Fig. 11 Comparison of time consumption betweenasynchronous policy update and synchronouspolicy update
图 12 考虑工况变化和不考虑工况变化统计结果对比
Fig. 12 Statistic results with and without consideration of dynamics of production condition
表 1 运行指标
Table 1 Operational indices
单元 运行指标 取值范围 (%) 竖炉 $a_1$: 磁管回收率 $a_{1\max} =84.8$ $a_{1\min} =81.3$ 磨矿单元1 $a_2$: 磨矿粒度 $a_{2\max} =84.0$ $a_{2\min} =48.6$ 磨矿单元2 $a_3$: 磨矿粒度 $a_{3\max} =88.8$ $a_{3\min} =63.3$ 强磁选 $a_4$: 精矿品位 $a_{4\max} =53.4$ $a_{4\min} =45.9$ $a_5$: 尾矿品位 $a_{5\max} =23.2$ $a_{5\min} =17.9$ 弱磁选 $a_6$: 精矿品位 $a_{6\max} =57.8$ $a_{6\min} =53.5$ $a_7$: 尾矿品位 $a_{7\max} =20.2$ $a_{7\min} =15.9$ 
下载: 导出CSV
-
[1] 柴天佑. 生产制造全流程优化控制对控制与优化理论方法的挑战. 自动化学报, 2009, 35(6): 641-649 doi: 10.3724/SP.J.1004.2009.00641Chai Tian-You. Challenges of optimal control for plant-wide production processes in terms of control and optimization theories. Acta Automatica Sinica, 2009, 35(6): 641-649 doi: 10.3724/SP.J.1004.2009.00641 [2] 丁进良, 杨翠娥, 陈远东, 柴天佑. 复杂工业过程智能优化决策系统的现状与展望. 自动化学报, 2018, 44(11): 1931-1943Ding Jin-Liang, Yang Cui-E, Chen Yuan-Dong, Chai Tian-You. Research progress and prospects of intelligent optimization decision making in complex industrial process. Acta Automatica Sinica, 2018, 44(11): 1931-1943 [3] 柴天佑, 丁进良, 王宏, 苏春翌. 复杂工业过程运行的混合智能优化控制方法. 自动化学报, 2008, 34(5): 505−515Chai Tian-You, Ding Jin-Liang, Wang Hong, Su Chun-Yi. Hybrid intelligent optimal control method for operation of complex industrial processes. Acta Automatica Sinica, 2008, 34(5): 505−515 [4] Huang X, Chu Y, Hu Y, Chai T. Production process management system for production indices optimization of mineral processing. IFAC Proceedings Volumes, 2005, 38(1): 178−183 [5] Ochoa S, Wozny G, Repke J U. Plantwide optimizing control of a continuous bioethanol production process. Journal of process Control, 2010, 20(9): 983−998 doi: 10.1016/j.jprocont.2010.06.010 [6] Ding J, Chai T, Wang H, Wang J, Zheng X. An intelligent factory-wide optimal operation system for continuous production process. Enterprise Information Systems, 2016, 10(3): 286−302 doi: 10.1080/17517575.2015.1065346 [7] Ding J, Modares H, Chai T, Lewis F L. Data-based multiobjective plant-wide performance optimization of industrial processes under dynamic environments. IEEE Transactions on Industrial Informatics, 2016, 12(2): 454−465 doi: 10.1109/TII.2016.2516973 [8] Chai T, Ding J, Wang H. Multi-objective hybrid intelligent optimization of operational indices for industrial processes and application. IFAC Proceedings Volumes, 2011, 44(1): 10517−10522 doi: 10.3182/20110828-6-IT-1002.01753 [9] Ding J, Yang C, Chai T. Recent progress on data-based optimization for mineral processing plants. Engineering, 2017, 3(2): 183−187 doi: 10.1016/J.ENG.2017.02.015 [10] Li J, Ding J, Chai T, Lewis F L. Nonzero-sum game reinforcement learning for performance optimization in large-scale industrial processes. IEEE Transactions on Cybernetics, 2019, 50(9): 4132−4145 [11] Liu C, Ding J, Sun J. Reinforcement learning based decision making of operational indices in process industry under changing environment. IEEE Transactions on Industrial Informatics, 2021, 17(4): 2727−2736 doi: 10.1109/TII.2020.3005207 [12] Lewis F L, Vrabie D, Vamvoudakis K. Reinforcement learning and feedback control. IEEE Control Systems, 2012, 32(6): 76−105 doi: 10.1109/MCS.2012.2214134 [13] Bertsekas D P, Tsitsiklis J N. Neuro-Dynamic Programming. Nashua: Athena Scientific, 1996. [14] Bertsekas D P. Proper policies in infinite-state stochastic shortest path problems. IEEE Transactions on Automatic Control, 2018, 63(11): 3787−3792 doi: 10.1109/TAC.2018.2811781 [15] Liu D, Wang D, Li H. Decentralized stabilization for a class of continuous-time nonlinear interconnected systems using online learning optimal control approach. IEEE Transactions on Neural Networks and Learning Systems, 2013, 25(2): 418−428 [16] Na J, hao J, Gao G, Li Z. Output-feedback robust control of uncertain systems via online data-Driven learning. IEEE Transactions on Neural Networks and Learning Systems, 2020, 32(6): 2650−2662 [17] Song R, Lewis F L, Wei Q. Off-policy integral reinforcement learning method to solve nonlinear continuous-time multiplayer nonzero-sum games. IEEE Transactions on Neural Networks and Learning Systems, 2016, 28(3): 704−713 [18] Modares H, Nageshrao S P, Lopes G A D, Babuska R, Lewis F L. Optimal model-free output synchronization of heterogeneous systems using off-policy reinforcement learning. Automatica, 2016, 71: 334−341 doi: 10.1016/j.automatica.2016.05.017 [19] Bertsekas D P. Multiagent reinforcement learning: rollout and policy iteration. IEEE/CAA Journal of Automatica Sinica, 2021, 8(2): 249−272 doi: 10.1109/JAS.2021.1003814 [20] Liang M, Wang D, Liu D. Neuro-optimal control for discrete stochastic processes via a novel policy iteration algorithm. IEEE Transactions on Systems, Man, and Cybernetics: Systems, 2019, 50(11): 3972−3985 [21] Zhang H, Luo Y, Liu D. Neural-network-based near-optimal control for a class of discrete-time affine nonlinear systems with control constraints. IEEE Transactions on Neural Networks, 2009, 20(9): 1490−1503 doi: 10.1109/TNN.2009.2027233 [22] Marvi Z, Kiumarsi B. Safe reinforcement learning: a control barrier function optimization approach. International Journal of Robust and Nonlinear Control, 2021, 31(6): 1923−1940 doi: 10.1002/rnc.5132 [23] Greene M L, Deptula P, Nivison S, Dixon W E. Sparse learning-based approximate dynamic programming with barrier constraints. IEEE Control Systems Letters, 2020, 4(3): 743−748 doi: 10.1109/LCSYS.2020.2977927 [24] Bellman R, Åström K J. On structural identifiability. Mathematical Biosciences, 1970, 7(3-4): 329−339 doi: 10.1016/0025-5564(70)90132-X [25] Luo B, Yang Y, Liu D. Policy iteration Q-learning for data-based two-player zero-sum game of linear discrete-time systems. IEEE Transactions on Cybernetics, 2021, 51(7): 3630−3640 doi: 10.1109/TCYB.2020.2970969 [26] Kiumarsi B, Lewis F L. Actor-critic-based optimal tracking for partially unknown nonlinear discrete-time systems. IEEE Transactions on Neural Networks and Learning Systems, 2014, 26(1): 140−151 [27] Zhang R, Tao J. Data-driven modeling using improved multi-objective optimization based neural network for coke furnace system. IEEE Transactions on Industrial Electronics, 2017, 64(4): 3147−3155 doi: 10.1109/TIE.2016.2645498 [28] Wang D, Ha M, Qiao J. Self-learning optimal regulation for discrete-time nonlinear systems under event-driven formulation. IEEE Transactions on Automatic Control, 2020, 65(3): 1272−1279 doi: 10.1109/TAC.2019.2926167 [29] Lewis F L, Liu D. Reinforcement Learning and Approximate Dynamic Programming for Feedback Control. New York: John Wiley & Sons, 2013. [30] Li J, Ding J, Chai T, Lewis F L, Jagannathan S. Adaptive interleaved reinforcement learning: robust stability of affine nonlinear systems with unknown uncertainty. IEEE Transactions on Neural Networks and Learning Systems, 2022, 33(1): 270-280 doi: 10.1109/TNNLS.2020.3027653 [31] 袁兆麟, 何润姿, 姚超, 李佳, 班晓娟. 基于强化学习的浓密机底流浓度在线控制算法. 自动化学报, 2021, 47(7): 1558-1571Yuan Zhao-Lin, He Run-Zi, Yao Chao, Li Jia, Ban Xiao-Juan. Online reinforcement learning control algorithm for concentration of thickener underflow. Acta Automatica Sinica, 2021, 47(7): 1558-1571 [32] Lowe R, Wu Y, Tamar A, Harb J, Abbeel P, Mordatch I. Multi-agent actor-critic for mixed cooperative-competitive environments. Advances in Neural Information Processing Systems, 2017, 6379-6390 [33] Sutton R S, Barto A G. Reinforcement Learning: An Introduction. Cambridge: MIT press, 2018. 期刊类型引用(1)
1. 李进，岳华峰，程生博，彭一帆，黄备备. 供应链质量追溯的烟草叶片图像帧特征动态识别方法. 计算技术与自动化. 2024(01): 111-116 . 
百度学术其他类型引用(6)
-
下载:
计量
- 文章访问数: 1464
- HTML全文浏览量: 203
- PDF下载量: 330
- 被引次数: 7
