A Modeling Method of Blast Furnace Burden Surface Accumulation Shape Based on the Motion Trajectory and Radial Distance
-
摘要: 高炉料面形貌是反映煤气流分布和煤气利用率的关键指标, 研究高炉料面炉料堆积形状数学建模方法对实现高炉精准布料控制和“双碳”战略在钢铁行业落地具有重要意义. 针对高炉多环布料情况下料面堆积形状预测难的问题, 本文提出了一种基于炉料运动轨迹和径向移动距离的高炉料面炉料堆积形状建模方法. 首先, 提出了一种与炉料初始状态和溜槽状态相关的炉料运动轨迹建模方法, 获取炉料从节流阀至料面的炉料运动轨迹, 并确定炉料在炉喉空区的内轨迹曲线和外轨迹曲线. 然后, 基于炉料运动轨迹和初始料面形状, 以体积守恒原则为约束, 提出了一种基于炉料径向移动距离的高炉料面炉料堆积形状数学建模方法, 获取炉料在料面的堆积形状. 最后, 基于某钢铁厂2# 高炉的尺寸建立离散单元法 (Discrete element method, DEM) 仿真模型, 模型仿真结果验证了所提方法的准确性和有效性.Abstract: The blast furnace burden surface is the key index to reflect the distribution of gas flow and the utilization rate of gas. Studying the mathematical modeling method of burden flow accumulation shape on the blast furnace burden surface is of great significance to realize the precise charging control and the implementation of “dual carbon” strategy in the steel industry. Aiming at the difficulty of predicting the burden flow accumulation shape in the blast furnace multi-ring charging, a modeling method for the accumulation shape of the burden flow on the blast furnace burden surface based on the burden flow motion trajectory and radial movement distance is proposed. Firstly, a modeling method of burden flow motion trajectory relate to the burden flow state and chute state is proposed to obtain the motion trajectory of burden flow from throttle valve to the burden surface, and further determine the internal and external trajectory of burden flow in the blast throat. Secondly, a mathematical modeling method of burden flow accumulation on the blast furnace burden surface based on the radial moving distance is proposed to obtain the accumulation shape of burden flow on the burden surface according to the motion trajectory, initial shape of burden surface, and the principle of volume conservation. Finally, a discrete element method (DEM) simulation model is established based on the 2# blast furnace of a steel plant, and the simulation results verify the accuracy and effectiveness of the proposed method.
-
现代控制系统变得越来越复杂, 相应地, 控制系统的安全性和可靠性就显得尤为重要. 但是, 实际系统运行时发生的各种故障会降低其可靠性甚至破坏系统的稳定性从而造成严重的安全事故. 因此, 为了提高控制系统的可靠性和安全性, 故障诊断技术在过去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 Schematic diagram of the moving process of burden flow on the blast furnace top
图 3 炉料与溜槽碰撞前后速度关系示意图
Fig. 3 Schematic diagram of the velocity relationship between the burden flow and chute collision
图 4 炉料在溜槽上位置示意图
Fig. 4 Schematic diagram of the position of the burden flow on the chute
图 5 炉料在料面落点位置和速度分布示意图
Fig. 5 Schematic diagram of the position and velocity distribution of the burden flow on the burden surface
图 9 传感器安装位置与料面堆积截面
Fig. 9 Installation positions of sensors and accumulation section of burden surface
图 10 基于DEM仿真的高炉料面堆积形状
Fig. 10 The accumulation shape of blast furnace burden surface based on DEM simulation
图 11 多环布料料面形状堆积更新过程
Fig. 11 The update process of burden surface accumulation shape of the multi-loop charging
图 12 基于DEM仿真和数学模型预测的料面堆积形状
Fig. 12 The burden surface accumulation shape based on DEM simulation and mathematical model
图 13 基于DEM仿真和数学模型预测的料面堆积形状绝对误差
Fig. 13 The absolute error of burden surface accumulation shape based on DEM simulation and mathematical model
表 1 DEM仿真中的粒子属性
Table 1 The particle properties used in DEM simulation
参数 数值 焦炭属性 半径(m) 0.04 ~ 0.06 泊松比(P−P) 0.22 剪切模量(Pa) 2.2×107 密度(kg/m3) 1050 恢复系数 碰撞恢复系数(P−P; P−C) 0.18; 0.2 静摩擦系数(P−P; P−C) 0.56; 0.41 滚动摩擦系数(P−P; P−C) 0.15; 0.09 (P−P和P−C分别表示粒子与粒子、粒子与溜槽接触) 
下载: 导出CSV
表 2 仿真中高炉布料操作参数
Table 2 The charging operation parameters in the simulation
布料环 环1 环2 环3 环4 环5 溜槽转速(°/s) 60 60 60 60 60 溜槽倾角(°) 41 39 36 33 30 旋转圈数 4 3 3 2 1
下载: 导出CSV
表 3 基于数学模型的料面高度预测性能
Table 3 The charging operation parameters in the burden furnace simulation
布料环 环1 环2 环3 环4 环5 合计轮廓点数 料面轮廓点 14 19 25 31 34 123 RMSE 0.0429 0.0240 0.0214 0.0259 0.0376 0.0308 $\Delta h \le 0.08$m $H(n)$ 14 19 25 31 32 121 ${\rm{HR}}$ 100% 100% 100% 100% 94.12% 98.37% $\Delta h \le 0.06$m $H(n)$ 11 19 25 31 32 116 ${\rm{HR}}$ 78.57% 100% 100% 96.77% 91.18% 94.31%
下载: 导出CSV
-
[1] 蒋珂, 蒋朝辉, 谢永芳, 潘冬, 桂卫华. 基于动态注意力深度迁移网络的高炉铁水硅含量在线预测方法. 自动化学报, 2023, 49(5): 949−963Jiang Ke, Jiang Zhao-Hui, Xie Yong-Fang, Pan Dong, Gui Wei-Hua. Online prediction method for silicon content of molten iron in blast furnace based on dynamic attention deep transfer network. Acta Automatica Sinica, 2023, 49(5): 949−963 [2] 周平, 张丽, 李温鹏, 戴鹏, 柴天佑. 集成自编码与PCA的高炉多元铁水质量随机权神经网络建模. 自动化学报, 2018, 44(10): 1799-1811Zhou Ping, Zhang Li, Li Wen-Peng, Dai Peng, Chai Tian-You. Modeling of blast furnace multi-element molten iron quality with random weight neural network based on self-encoding and PCA. Acta Automatica Sinica, 2018, 44(10): 1799-1811 [3] Jimenez J, Mochon J, Formoso A, de Ayala J S. Burden distribution analysis by digital image processing in a scale model of a blast furnace shaft. ISIJ International, 2000, 40(2): 114-120 doi: 10.2355/isijinternational.40.114 [4] Mitra T, Saxén H. Simulation of burden distribution and charging in an ironmaking blast furnace. IFAC-PapersOnLine, 2015, 48(17): 183-188 doi: 10.1016/j.ifacol.2015.10.100 [5] Kajiwara Y, Jimbo T, Joko T, Aminaga Y, Inada T. Investigation of bell-less charging based on full scale model experiments. Transactions of the Iron and Steel Institute of Japan, 1984, 24(10): 799-807 doi: 10.2355/isijinternational1966.24.799 [6] Mio H, Komatsuki S, Akashi M, Shimosaka A, Shirakawa Y, Hidaka J, Kadowaki M, Matsuzaki S, Kunitomo K. Effect of Chute Angle on Charging Behavior of Sintered Ore Particles at Bell-less Type Charging System of Blast Furnace by Discrete Element Method. ISIJ International, 2009, 49(4): 479-486 doi: 10.2355/isijinternational.49.479 [7] Kou M Y, Wu S L, Zhou H, Yu Y M, Xu J. Numerical Investigation of Coke Collapse and Size Segregation in the Bell-less Top Blast Furnace. Isij International, 2018, 58(11): 2018-2024 doi: 10.2355/isijinternational.ISIJINT-2018-415 [8] Zhou K, Jiang Z H, Pan D, Gui W H, Huang J C. Influence of charging parameters on the burden flow velocity and distribution on the blast furnace chute based on discrete element method. Steel Research International, 2022, 93(1): Article No. 2100332 [9] Hong Z B, Zhou H, Wu J L, Zhan L L, Fan Y B, Zhang Z W, et al. Effects of operational parameters on particle movement and distribution at the top of a bellless blast furnace based on discrete element method. Steel Research International, 2021, 92(1): Article No. 2000262 [10] 马洪佑, 王振阳, 戴建华, 袁军, 李秀亮, 王永龙. 不同溜槽形状下的料流偏析现象. 钢铁, 2020, 55(09): 23-28Ma Hong-You, Wang Zhen-Yang, Dai Jian-Hua, Yuan Jun, Li Xiu-Liang, Wang Yong-Long.Segregation of material flow under different chute shapes.Iron and Steel, 2022, 55(09): 23-28 [11] 孙俊杰, 狄瞻霞, 李家新, 卢开成. 溜槽形状及倾角对料流运动的影响. 钢铁, 2019, 54(04): 19-23Sun Jun-Jie, Di Zhan-Xia, Li Jia-Xin, Lu Kai-Cheng. Influence of inclination and shape of chute on movement of burden flow.Iron and Steel, 2019, 54(04): 19-23 [12] Kou M Y, Xu J, Wu S L, Zhou H, Gu K, Yao S, Wen B J. Effect of cross-section shape of rotating chute on particle movement and distribution at the throat of a bell-less top blast furnace. Particuology, 2019, 44: 194-206 doi: 10.1016/j.partic.2018.07.010 [13] Govender N, Wilke D N, Wu C Y, Tuzun U, Kureck H. A numerical investigation into the effect of angular particle shape on blast furnace burden topography and percolation using a GPU solved discrete element model. Chemical Engineering Science, 2019, 204: 9-26 doi: 10.1016/j.ces.2019.03.077 [14] 寇明银, 吴胜利, 周恒, 于亚光, 顾凯. 高炉DEM模拟用炉料系数的测定及其应用. 钢铁, 2018, 53(12): 30-36+61 doi: 10.13228/j.boyuan.issn0449-749x.20180163Kou Ming-Yin, Wu Sheng-Li, Zhou Heng, Yu Ya-Guang, Gu Kai.Measurements and application of burden coefficients for DEM simulation in blast furnace. Iron and Steel, 2018, 53(12): 30-36+61 doi: 10.13228/j.boyuan.issn0449-749x.20180163 [15] 夏修浩, 周连勇, 马华庆, 赵永志. 颗粒形状模型对高炉布料过程DEM模拟的影响. 钢铁研究学报, 2021, 33(12): 1228-1236Xia Xiu-Hao, Zhou Lian-Yong, Ma Hua-Qing, Zhao Yong-Zhi. Effect of particle shape model on DEM simulation of charging process in blast furnace.Journal of Iron and Steel Research, 2021, 33(12): 1228-1236 [16] Mio H, Narita Y, Matsuzaki S, Nishioka K, Nomura S. Measurement of particle charging trajectory via rotating chute of 1/3-scale blast furnace and its comparing with numerical analysis using Discrete Element Method. Powder Technology, 2019, 344: 797-803 doi: 10.1016/j.powtec.2018.12.047 [17] Wei S Y, Wei H, Saxen H, Yu Y W. Numerical analysis of the relationship between friction coefficient and repose angle of blast furnace raw materials by discrete element method. Materials, 2022, 15(3): Article No. 903 [18] Holzinger G, Schatzl M. Effect of chute start angle and hopper change on burden distribution during the charging process of a bell-less top blast furnace with two parallel hoppers. Powder Technology, 2022, 395: 669-680 doi: 10.1016/j.powtec.2021.10.005 [19] Yu Y W, Saxen H. Particle Flow and Behavior at Bell-Less Charging of the Blast Furnace. Steel Research International, 2013, 84(10): 1018-1033 [20] Mitra T, Saxén H. Discrete element simulation of charging and mixed layer formation in the ironmaking blast furnace. Computational Particle Mechanics, 2016, 3(4): 541-555 doi: 10.1007/s40571-015-0084-1 [21] Mitra T, Saxen H. Investigation of Coke Collapse in the Blast Furnace Using Mathematical Modeling and Small Scale Experiments. ISIJ International, 2016, 56(9): 1570-1579 doi: 10.2355/isijinternational.ISIJINT-2016-114 [22] Radhakrishnan V R, Ram K M. Mathematical model for predictive control of the bell-less top charging system of a blast furnace. Journal of Process Control, 2001, 11(5): 565-586 doi: 10.1016/S0959-1524(00)00026-3 [23] 朱清天, 程树森, 魏志江, 郭喜斌. 高炉炉料落点的确定. 中国冶金, 2006, (09): 24-26Zhu Qing-Tian, Cheng Shu-Sen, Wei Zhi-Jiang, Guo Xi-Bin. Establishment of Landing Point of Burden in Blast Furnace. China Metallurgy, 2006, (09): 24-26 [24] 杜鹏宇, 程树森, 胡祖瑞, 吴桐. 高炉无钟炉顶布料料流宽度数学模型及试验研究. 钢铁, 2010, 45(01): 14-18Du Peng-Yu, Cheng Shu-Sen, Hu Zu-Rui, Wu Tong. Mathematical Model of Burden Width in a Bell-Less Top Blast Furnace and Modeling Experimental Research. Iron and Steel, 2010, 45(01): 14-18 [25] Fu D, Chen Y, Zhou C Q. Mathematical modeling of blast furnace burden distribution with non-uniform descending speed. Applied Mathematical Modelling, 2015, 39(23-24): 7554-7567 doi: 10.1016/j.apm.2015.02.054 [26] 张森, 李酉, 陈先中, 尹怡欣. 高炉料面形状双驱动模型研究. 控制理论与应用, 2020, 37(05): 978-986Zhang Sen, Li You, Chen Xian-Zhong, Yin Yi-Xin. Research on double-driven model of blast furnace burden profile. Control Theory & Applications, 2020, 37(05): 978-986 [27] Fojtik D, Tuma J, Faruzel P. Computer modelling of burden distribution in the blast furnace equipped by a bell-less top charging system. Ironmaking & Steelmaking, 2021, 48(10): 1226-1238 [28] Nag S, Gupta A, Paul S, Gavel D J, Aich B. Prediction of Heap Shape in Blast Furnace Burden Distribution. ISIJ International, 2014, 54(7): 1517-1520 doi: 10.2355/isijinternational.54.1517 [29] Li M, Wei H, Ge Y, Xiao G C, Yu Y W. A Mathematical model combined with radar data for bell-less charging of a blast furnace. Processes, 2020, 8(2): Article No. 239 [30] 周科, 蒋朝辉, 桂卫华, 黄建才, 朱既承. 基于坐标变换的高炉U型溜槽上炉料运动轨迹数学建模研究. 2021 中国自动化大会. 昆明, 中国: 2021. 289−294Zhou Ke, Jiang Zhao-Hui, Gui Wei-Hua, Huang Jian-Cai, Zhu Ji-Cheng. Research on mathematical modeling of three-dimensional movement of particle on U-shaped chute of blast furnace based on coordinate transformation. In: Proceedings of the 2021 China Automation Conference Proceedings. Kunming, China: 2021. 289−294 [31] 冯建强, 孙诗一. 四阶龙格—库塔法的原理及其应用. 数学学习与研究, 2017, (17) : 3-5Feng Jian-Qiang, Sun Shi-Yi. Principle and application of fourth-order Runge-Kutta method. Mathematics Learning and Research, 2017, (17): 3-5 [32] 李晓理, 刘德馨, 周翔, 陈先中. 高炉布料设定值优化控制. 控制理论与应用, 2015, 32(12) : 1660-1668 doi: 10.7641/CTA.2015.40856Li Xiao-Li, Liu De-Xin, Zhou Xiang, Chen Xian-Zhong. Setting value optimal control for blast furnace burden distribution. Control Theory & Applications, 2015, 32(12): 1660-1668 doi: 10.7641/CTA.2015.40856 [33] 马财生, 任廷志. 无钟高炉炉料分布预测模型. 工程科学学报, 2017, 39(02): 276-282Ma Cai-Sheng, Ren Yan-Ting. Burden distribution prediction model in a blast furnace with bell-less top. Chinese Journal of Engineering, 2017, 39(02): 276-282 [34] Cundall P A, Strack O D L. A discrete numerical model for granular assemblies. Geotechnique, 1979, 29(1): 47-65 doi: 10.1680/geot.1979.29.1.47 [35] Yu Y W, Saxen H. Segregation behavior of particles in a top hopper of a blast furnace. Powder Technology, 2014, 262: 233-241 doi: 10.1016/j.powtec.2014.04.010 [36] Zhou K, Jiang Z, Pan D, Gui W H, Huang J C, Xu C. Research on the velocity distribution law of the coke in the chute of blast furnace based on discrete element method. Computational Particle Mechanics, 2022: 1-9 [37] Ester M, Kriegel H P, Sander J, et al. A density-based algorithm for discovering clusters in large spatial databases with noise. KDD, 1996: 226-231 期刊类型引用(1)
1. 李进，岳华峰，程生博，彭一帆，黄备备. 供应链质量追溯的烟草叶片图像帧特征动态识别方法. 计算技术与自动化. 2024(01): 111-116 . 
百度学术其他类型引用(6)
-
下载:
计量
- 文章访问数: 688
- HTML全文浏览量: 229
- PDF下载量: 206
- 被引次数: 7
