Riemannian Metric Based Performance Monitoring and Diagnosis for a Class of Feedback Control Systems
-
摘要: 针对复杂工业系统对性能衰退的容忍度低等问题, 提出基于系统性能预测的一类反馈控制系统过程监测方法, 通过黎曼度量对控制性能衰退程度进行预测与监测, 并给出发生故障的类型, 以提升过程监测系统的实时性与准确性. 首先, 利用系统的实时数据, 计算系统性能衰退的预测指标; 其次, 利用黎曼度量对系统性能衰退程度进行预测与监测, 并利用随机算法给出对应的阈值来诊断系统性能衰退; 最后, 通过计算各类引发系统性能衰退的故障的性能预测指标集合的中心和阈值, 实现故障的实时定位. 所提出的方法通过三容水箱仿真实验平台进行验证.Abstract: Associated with the issue that the industrial processes have low tolerance in performance degradation, a performance monitoring and diagnosis scheme for a class of feedback control systems is developed based on system performance prediction. The Riemannian metric is then adopted to predict and monitor the degree of the control performance degradation, and classify the fault type, so as to improve the real-time ability and accuracy of the performance monitoring systems. Firstly, a prediction indicator is introduced to predict the system performance degradation using online process data. Then, the Riemannian metric is implemented to predict and monitor the degree of the control performance degradation, and the randomized algorithm is adopted to set the threshold for diagnosing the performance degradation. Finally, by computing the center and threshold for the performance predictor set caused by all the possible faults, the fault localization can be achieved. The effectiveness of the proposed method is verified by a case study on a three-tank system.
-
Key words:
- Performance monitoring /
- Riemannian metric /
- performance diagnosis /
- fault location
-
随着现代工业复杂化程度的不断提升, 各级生产环节的关联也越来越密切, 导致系统变量间相互关联、相互耦合, 一旦某个环节发生性能衰退就可能随着链式反应在系统中传播, 进而导致整个系统瘫痪, 造成不可挽回的损失. 为保障产品生产的质量和生产效率, 提高工业生产过程的安全性和可靠性, 需要对工业生产过程或运行设备进行实时监测. 因此过程监测与故障诊断在研究和工业应用领域都受到了极大的关注[1-6], 成为自动控制领域中的一个研究热点, 对保障工业过程安全性、可靠性具有重要意义.
由于现代工业过程对系统性能、效率和可靠性要求的不断提高, 控制性能监测技术也得到了广泛关注[7-8]. 系统控制性能监测的核心思想是将系统实际运行状况与根据系统设计的基准进行对比, 判定当前系统是否运行良好. 目前为止研究的控制性能监测方法主要集中在对控制器的性能进行评价, 指出其与最优控制性能之间的差距[7-10]. 近年来, 涌现了不少控制性能监测的改进算法. 具体而言, 基于关键性能指标预测的控制性能监测方法得到了广泛关注[11-13]; Li等[14]提出了基于稳定性能衰退评估的故障诊断方法; Tao等[15]建立一个集成性能评估、故障检测和诊断的框架, 最大限度地利用现有监测数据来解决性能评估问题; Li等[16]借助数据驱动技术提出了基于性能衰退预测的反馈控制系统的性能监测与故障诊断方法. 一般而言, 控制性能指标可分为确定性指标、鲁棒性指标、随机性指标. 确定性指标包含了调节时间、衰减率等传统意义上的控制回路动态品质指标, 鲁棒性指标衡量模型失配下的系统稳定性和控制品质, 随机性指标则是一种统计意义下的指标, 可以用来衡量系统的动态变化. 在随机性能评价领域, 专家学者们做出了巨大贡献, 提出了基于历史数据的基准、最小方差基准、广义最小方差等准则的性能监测方法[17-20]. 目前大多数控制性能监测方法主要集中于系统部件故障的诊断, 而很少关注基于系统性能衰退程度的诊断. 这些性能衰退不仅可能由系统部件故障引起, 也可能由系统控制参数不匹配引起. 另一方面, 目前的控制性能监测的方法很少对系统的运行性能进行预测. 而基于系统性能预测的诊断可以及时地发现可能发生的问题, 进而为系统性能修复提供时间. 因此, 本文研究基于性能衰退预测的性能监测方法, 旨在实现实时性能诊断.
在故障诊断领域, 距离度量可有效地根据待测数据集和无故障运行数据集之间的相似程度来诊断故障的发生. 距离度量是一种空间属性, 选取不同的度量函数, 可以将数据集映射到不同的空间中, 直接关系到故障诊断结果. 因此, 度量函数的研究在故障诊断领域逐渐发展起来. 常见的距离度量有: 欧氏距离[21-23]、马氏距离[24-27]等. 近年来, 随着流形研究的飞速发展, 黎曼度量作为流形曲面上的测地线距离, 越来越受到人们的关注[28-29], 专家学者们也开始将黎曼度量应用于故障诊断中. An等[30]提出了一种基于黎曼度量和一维卷积神经网络的端到端无监督域自适应轴承故障诊断方法, 该方法具有较强的故障识别能力和领域不变性, 适用于频繁变化的工作环境. 周美含等[31]提出了一种基于黎曼度量的单基地雷达目标检测新方法, 通过计算噪声协方差矩阵的黎曼均值, 将其与接收信号协方差矩阵之间的黎曼度量作为检测统计量, 实现故障诊断. 与传统的基于欧氏度量的检测方法相比, 基于黎曼度量的方法显著提高了低信噪比和单快拍下的目标检测性能. 由于黎曼度量适用于变量间高度耦合的情况, 因此故障诊断领域内越来越多的专家学者开始采用黎曼度量作为工具度量数据间的距离[32-34]. 然而, 目前基于黎曼度量的故障诊断研究多针对静态系统数据展开, 而很少有将黎曼度量用于动态系统的性能监测与诊断中.
本文从性能预测角度出发, 针对一类带有反馈控制环节的动态系统, 提出一种基于黎曼度量的性能预测与评估指标, 并基于该指标实现系统性能监测与诊断. 由于在实际工业过程中, 系统的模型结构与内部参数往往存在着一些不确定的因素, 例如参数摄动、不确定性以及测量仪器精度造成的误差和模型线性化引起的不确定性等, 因此本文针对含不确定性的反馈控制系统展开研究.
本文结构如下: 第1节介绍具有模型不确定性的动态系统, 并引入黎曼度量, 提出基于黎曼度量的系统变化检测的基本思路; 第2节提出动态系统在反馈控制下的性能衰退预测指标, 并研究基于模型和数据驱动的性能预测指标的辨识方法, 然后给出基于黎曼度量的控制性能监测方法; 第3节基于黎曼度量提出对引发系统性能衰退的故障进行定位的方案; 第4节采用三容水箱系统对本文所提出的性能监测与故障定位方案进行验证; 最后总结本文. 与现有的诊断方法相比, 本文具有如下创新性:
1)提出了反馈控制系统的性能衰退预测方法并给出了对应的在线辨识方案;
2)基于黎曼度量实现了动态系统的控制性能实时监测与诊断. 与现有的基于黎曼度量的故障检测方法直接利用采集的测量数据形成正定矩阵不同, 本文通过在线辨识的性能预测矩阵来监测系统性能的衰退程度. 因此, 本文所提出的方法可用于动态系统的性能监测与诊断.
1. 系统介绍与问题描述
本节首先介绍论文中考虑的反馈控制系统与基于黎曼度量的系统变化检测方法, 在此基础上提出本文拟解决的问题.
1.1 反馈控制系统
考虑如下线性时不变系统
$$ {\boldsymbol x}(k + 1) = {\boldsymbol A}{\boldsymbol x}(k) + {\boldsymbol B}{\boldsymbol u}(k)+{\boldsymbol w}(k)\\ $$ (1) 其中, $ {\boldsymbol x}(k) \in {{\bf R}^n} $表示系统的状态变量, $ {\boldsymbol u}(k)\in {{\bf R}^m} $表示控制输入. 被控输出$ {\boldsymbol y}(k)\in {{\bf R}^p} $可描述为
$$ \begin{array}{l} {\boldsymbol y}(k) = {\boldsymbol C}{\boldsymbol x}(k) + {\boldsymbol D}{\boldsymbol u}(k)+{\boldsymbol\eta}(k) \end{array} $$ (2) 在式(1)和式(2)中, $ {\boldsymbol w}(k),{\boldsymbol\eta}(k) $分别表示零均值的系统过程与测量噪声, 满足
$$ \begin{array}{l} {\cal{E}}\left( \left[ {\begin{array}{*{20}{c}} {\boldsymbol w}(i) \\ {\boldsymbol\eta}(i) \\ {\boldsymbol x}(0) \end{array}} \right] \left[ {\begin{array}{*{20}{c}} {\boldsymbol w} (j) \\ {\boldsymbol \eta}(j)\\ {\boldsymbol x}(0) \end{array}} \right]^{\rm T} \right) = \left[ {\begin{array}{*{20}{c}} \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{\Sigma}}}_w & {\bf 0} \\ {\bf 0} & {\boldsymbol{{ \Sigma}}}_\eta \end{array}} \right]& {\bf 0} \\ {\bf 0} & {\boldsymbol{{ \Pi}}}_0 \end{array}} \right] \end{array} $$ (3) 在式(3)中, $ {\boldsymbol \Sigma}_w $, $ {\boldsymbol \Sigma}_\eta $, $ {\boldsymbol \Pi}_0 $为已知的正定矩阵. $ {{\boldsymbol A}} $, $ {{\boldsymbol B}} $, $ {{\boldsymbol C}} $, $ {{\boldsymbol D}} $表示随不确定性参数变化的系统矩阵, 即
$$\begin{split} &\left[ {\begin{array}{*{20}{c}} {\boldsymbol A} & {\boldsymbol B} \\ {\boldsymbol C} & {\boldsymbol D} \end{array}} \right]=\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol A}}_{0}} & {{{\boldsymbol B}}_{0}} \\ {{{\boldsymbol C}}_{0}} & {{{\boldsymbol D}}_{0}} \end{array}} \right]+{\boldsymbol \Delta } \\ & {\boldsymbol \Delta} =\sum_{i=1}^{l} \delta_i \left[ {\begin{array}{*{20}{c}} {{\boldsymbol \Delta}_{{A_i}}} & {{\boldsymbol \Delta}_{{B_i}}} \\ {{\boldsymbol \Delta}_{{C_i}}} & {{\boldsymbol \Delta}_{{D_i}}} \end{array}} \right] \end{split}$$ (4) 其中, $ {{\boldsymbol A}_0} $, $ {{\boldsymbol B}_0} $, $ {{\boldsymbol C}_0} $, $ {{\boldsymbol D}_0} $表示已知的标称系统矩阵, $ {\boldsymbol \Delta} $表示不确定性, $ {\boldsymbol \Delta}_{{A_i}} $, $ {\boldsymbol \Delta}_{{B_i}} $, $ {\boldsymbol \Delta}_{{C_i}} $, $ {\boldsymbol \Delta}_{{D_i}} $为已知矩阵. $ \delta_i,\;i = 1,\cdots ,\;l $为在区间 $ [\bar{\delta}_{i,1},\bar{\delta}_{i,2}] $上均匀分布的随机变量.
本文考虑如下反馈控制
$$ \begin{array}{l} {\boldsymbol u}(k) = {\boldsymbol K}{\boldsymbol x}(k)+{\boldsymbol v}(k) \end{array} $$ (5) 其中, $ {\boldsymbol K} $为反馈控制器, $ {\boldsymbol v}(k)\in {{\bf R}^m} $为可由式(6)和式(7)表示的参考输入
$$ \begin{array}{l} {\boldsymbol x}_v(k+1) = {\boldsymbol A}_v{\boldsymbol x}_v(k)+{\boldsymbol B}_v{\boldsymbol r}(k) \end{array} $$ (6) $$ \begin{array}{l} {\boldsymbol v}(k) = {\boldsymbol C}_v{\boldsymbol x}_v(k) \quad\qquad\qquad\qquad\qquad \end{array} $$ (7) 其中, $ {\boldsymbol r}(k)\in {\bf R}^{k_r} $表示常数参考信号, $ {\boldsymbol x}_v(k)\in {\bf R}^{k_v} $表示对应的状态, $ {\boldsymbol A}_v,{\boldsymbol B}_v,{\boldsymbol C}_v $表示相应的矩阵.
1.2 基于黎曼度量的变化检测
本节介绍黎曼度量和基于黎曼度量的变化检测方法.
定义$ {\cal{P}}(m) $为所有$ m\times m $维的对称正定矩阵$ {\boldsymbol P} $所组成的集合. 该集合形成了${m(m+1)}/{2}$维的黎曼流形. 黎曼度量表示两点在流形曲面上的最短距离(测地线距离)[28-29, 35]. 给定集合$ {\cal{P}} $中的任意两个矩阵$ {\boldsymbol P}_1,{\boldsymbol P}_2 $, 其黎曼度量定义为
$$ \begin{split} d_R({\boldsymbol P}_1,{\boldsymbol P}_2) = \;& \left\Vert \ln {\boldsymbol P}_1^{-\frac{1}{2}}{\boldsymbol P}_2 {\boldsymbol P}_1^{-\frac{1}{2}}\right\Vert_{\rm{F}} =\\& \sqrt{ \sum_{j = 1}^m \left(\ln \lambda_j \left({\boldsymbol P}_1^{-1}{\boldsymbol P}_2 \right)\right)^2} \end{split} $$ (8) 其中, $\lambda_j \left({\boldsymbol P}_1^{-1}{\boldsymbol P}_2\right)$为${\boldsymbol P}_{1}^{-1}{{{\boldsymbol P}}_{2}}$的第j个特征值, $j = 1,\cdots ,\;m$. $\left\Vert \cdot \right\Vert_{\rm{F}}$表示Frobenious范数.
假设矩阵${\boldsymbol P}_{1},{\boldsymbol P}_2$是$m(m+1)/{2}$维流形上的点, 而黎曼度量可以有效地刻画两个正定矩阵之间的差异. 因此, 正定矩阵的黎曼度量可作为监测量用于检测系统可能发生的变化. 假设在理想或系统未发生变化的状态下, 过程数据集${ \Omega}_i, i = 1,\cdots ,\;M_0$被记录下来, 其中, $M_0$为数据集的个数. 提取数据集${ \Omega}_i$对应的正定特征矩阵${\boldsymbol P}_i$, 则可计算不同数据集间的黎曼度量.
在介绍基于黎曼度量的系统变化检测方法之前, 首先引入如下黎曼均值的概念[36-38]
$$ {\boldsymbol P}_z = \arg \min\limits_{{\boldsymbol P}_j, j = 1,\cdots ,\; M_0} \sum\limits_{i = 1}^{M_0} d_R^2({\boldsymbol P}_i,{\boldsymbol P}_j) $$ (9) 即离流形上所有点$ {\boldsymbol P}_i,i = 1,\cdots ,\;M_0 $的黎曼度量和最小的点就表示该流形的黎曼均值.
本文中将流形上任一点$ {\boldsymbol P}_i $与黎曼均值$ {\boldsymbol P}_z $之间的黎曼度量$ d_R^2({\boldsymbol P}_i,{\boldsymbol P}_z) $作为监测量用于变化检测, 并设定如下阈值
$$ \begin{array}{l} J_{{\rm{th}}} = \max\limits_{i\in \{1,\cdots ,\;M_0\}} d_R^2({\boldsymbol P}_i,{\boldsymbol P}_z) \end{array} $$ (10) 当采集到实时的系统数据后, 可提取新数据集的特征矩阵$ {\boldsymbol P}_{{\rm{new}}} $. 考虑到$ {\boldsymbol P}_{{\rm{new}}} $可以表征系统发生的变化, 因此可借助黎曼度量
$$ \begin{array}{l} J = d_R^2({\boldsymbol P}_{{\rm{new}}},{\boldsymbol P}_z) \end{array} $$ (11) 通过如下决策逻辑实现系统变化检测(如图1所示)
$$ \left\{ \begin{aligned} &J_{{\rm{th}}}-J\ge 0 \Longrightarrow {\text{无变化}}\\ &J_{{\rm{th}}}-J< 0 \Longrightarrow {\text{系统发生变化}} \end{aligned}\right. $$ (12) 1.3 问题描述
从式(8)和式(11)可知, 基于黎曼度量的变化检测方法可以有效地监测正定矩阵在方向和幅值方面的变化[38]. 对于反馈控制系统而言, 这些变化体现在由参数变化或乘性故障引起的系统动态特性的改变. 目前的过程监测方法多采用二次型性能指标, 这些性能指标为标量, 对参数变化或乘性故障引起的系统动态变化不够灵敏. 因此, 本文采用黎曼度量来实现反馈控制系统的性能监测.
一般而言, 反馈控制律的设计用于保证系统具有理想的性能(如能耗、时间、控制性能等). 不失一般性, 本文引入如下的性能指标来衡量系统的控制性能
$$ V(k) = \sum\limits_{i = k}^\infty \gamma^{i-k} {\cal{E}}\left({{{\boldsymbol y}^{\rm T}}(i){\boldsymbol Q}{\boldsymbol y}(i) + {{\boldsymbol u}^{\rm T}}(i){\boldsymbol R}{\boldsymbol u}(i)}\right) $$ (13) 其中, $ {\boldsymbol Q} \ge 0,{\boldsymbol R}> 0 $为加权矩阵, 可根据系统对性能的要求进行选取. $ {\cal{E}} $表示信号的均值. $ 0<\gamma< 1 $表示衰减因子.
很显然对于给定的控制律, 性能指标(13)都对应一个给定的值. 从应用的角度而言, 我们通常要保证系统的运行满足需求的性能. 然而实际中, 系统操作点或参数的改变、系统部件的老化、控制参数的不匹配等都可能引发系统性能衰退, 从而改变性能指标(13)的值. 因此, 需要设计合适的方法来预测系统性能的变化并将它应用于控制性能监测. 本文的目标是研究基于黎曼度量的系统控制性能监测与诊断方法. 具体而言, 本文主要进行以下工作:
1)选择合适的性能指标函数来预测系统性能;
2)基于性能预测指标, 运用黎曼度量对系统性能进行在线监测;
3)建立故障模态数据库, 对所发生系统性能衰退的类型进行识别.
2. 控制性能衰退预测与监测
基于黎曼度量的控制监测方法包括系统控制性能衰退预测、基于随机算法的阈值设定以及基于黎曼度量的控制性能监测.
2.1 控制性能衰退预测方法
为了实现控制性能衰退预测, 首先证明如下定理.
定理 1. 考虑系统(2), 给定反馈控制律(6), 则性能指标(13) 可表示成如下参数化形式
$$ V(k) = {\cal{E}} {\bar{{\boldsymbol x}}^{\rm T}(k)}{\boldsymbol P}{\bar{{\boldsymbol x}}(k)}+c $$ (14) 其中, $\bar{{\boldsymbol x}}(k)=\left[ {\boldsymbol x}^{\rm T}(k) \;\;{\boldsymbol x}_v^{\rm T}(k) \;\; {\boldsymbol r}^{\rm T}(k)\right]^{\rm T}$, $ {\boldsymbol P} $是正定矩阵, $ c $为正常数且满足
$$ {\boldsymbol P} = \bar{{\boldsymbol C}}_K^{\rm T}{\boldsymbol Q}\bar{{\boldsymbol C}}_K+\bar{{\boldsymbol K}}^{\rm T}{\boldsymbol R}\bar{{\boldsymbol K}}+\gamma \bar{{\boldsymbol A}}_K^{\rm T} {\boldsymbol P} \bar{{\boldsymbol A}}_K \quad\quad$$ (15) $$ c = \frac{{\rm{tr}}({\boldsymbol \Sigma}_\eta {\boldsymbol Q}) +\gamma {\rm{tr}}({\boldsymbol \Sigma}_w \bar{{\boldsymbol I}}{\boldsymbol P}\bar{{\boldsymbol I}}^{\rm T})}{1-\gamma }\;\; \quad\quad\quad\qquad $$ (16) 其中,
$$\begin{split} &\bar{{\boldsymbol A}}_K=\left[ {\begin{array}{*{20}{c}} {\boldsymbol A}+{\boldsymbol B}{\boldsymbol K} & {\boldsymbol B}{\boldsymbol C}_v & {\bf 0}\\ {\bf 0} & {\boldsymbol A}_v & {\boldsymbol B}_v\\ {\bf 0} & {\bf 0} & {\boldsymbol I} \end{array}} \right],\bar{{\boldsymbol I}}=\left[ {\begin{array}{*{20}{c}} {\boldsymbol I} \\ {\bf 0} \\ {\bf 0} \end{array}} \right]\notag\\ &\bar{{\boldsymbol C}}_K=\left[ {\begin{array}{*{20}{c}} {\boldsymbol C}+{\boldsymbol D}{\boldsymbol K} & {\boldsymbol D}{\boldsymbol C}_v &{\bf 0} \end{array}} \right], \bar{{\boldsymbol K}}=\left[ {\begin{array}{*{20}{c}} {\boldsymbol K}& {\boldsymbol C}_v & {\bf 0} \end{array}} \right]\notag \end{split}$$ 证明. 对于性能指标(13)而言, 可以得到如下递推方程, 即Bellman方程
$$ \begin{split} V(k) =\;& {\cal{E}}\left({\boldsymbol y}^{\rm T}(k){\boldsymbol Q}{\boldsymbol y}(k)+ {\boldsymbol u}^{\rm T}(k){\boldsymbol R}{\boldsymbol u}(k)\right)+\\ &\gamma V(k+1) \end{split} $$ (17) 考虑到闭环系统的状态空间可描述为
$$ \bar{{\boldsymbol x}}(k+1) = \bar{{\boldsymbol A}}_K\bar{{\boldsymbol x}}(k)+\bar{{\boldsymbol I}}{\boldsymbol w}(k) $$ (18) $$ {\boldsymbol y}(k) = \bar{{\boldsymbol C}}_K\bar{{\boldsymbol x}}(k)+{\boldsymbol\eta}(k)\quad\quad\;\; $$ (19) $$ {\boldsymbol u}(k) = \bar{{\boldsymbol K}} \bar{{\boldsymbol x}}(k)\qquad\qquad\qquad \; $$ (20) 则瞬时性能指标可表示为
$$ \begin{split} {\boldsymbol y}^{\rm T}(k)&{\boldsymbol Q}{\boldsymbol y}(k)+{\boldsymbol u}^{\rm T}(k){\boldsymbol R}{\boldsymbol u}(k)\notag =\\ &\bar{{\boldsymbol x}}^{\rm T}(k)\left(\bar{{\boldsymbol C}}_K^{\rm T}{\boldsymbol Q}\bar{{\boldsymbol C}}_K+\bar{{\boldsymbol K}}^{\rm T}{\boldsymbol R}\bar{{\boldsymbol K}}\right){\bar{{\boldsymbol x}}(k)}\notag\;+\\ &2\bar{{\boldsymbol x}}^{\rm T}(k)\bar{{\boldsymbol C}}_K^{\rm T}{\boldsymbol Q}{\boldsymbol\eta}(k)+{\boldsymbol\eta}^{\rm T}(k){\boldsymbol Q}{\boldsymbol\eta}(k) \end{split} $$ 将式(14)代入式(17), 可得
$$ \begin{split} & {\cal{E}} \left({\bar{{\boldsymbol x}}^{\rm T}(k)}{\boldsymbol P}{\bar{{\boldsymbol x}}(k)}\right)+c\notag=\\&\qquad {\cal{E}} \left({\boldsymbol y}^{\rm T}(k){\boldsymbol Q}{\boldsymbol y}(k)+{\boldsymbol u}^{\rm T}(k){\boldsymbol R}{\boldsymbol u}(k)\right)\notag+\\&\qquad\gamma {\cal{E}}\left( {\bar{{\boldsymbol x}}^{\rm T}(k+1)}{\boldsymbol P}{\bar{{\boldsymbol x}}(k+1)}+c\right)\notag=\\&\qquad{\cal{E}} \left( \bar{{\boldsymbol x}}^{\rm T}(k)\left(\bar{{\boldsymbol C}}_K^{\rm T}{\boldsymbol Q}\bar{{\boldsymbol C}}_K+\bar{{\boldsymbol K}}^{\rm T}{\boldsymbol R}\bar{{\boldsymbol K}}\right){\bar{{\boldsymbol x}}(k)}\right)+\\&\qquad{\cal{E}} \left( {\boldsymbol\eta}^{\rm T}(k){\boldsymbol Q}{\boldsymbol\eta}(k)\right)+\gamma {\cal{E}} \left( \bar{{\boldsymbol x}}^{\rm T}(k)\bar{{\boldsymbol A}}_K^{\rm T} {\boldsymbol P} \bar{{\boldsymbol A}}_K {\bar{{\boldsymbol x}}(k)}\right.+\\&\qquad2\bar{{\boldsymbol x}}^{\rm T}(k)\bar{{\boldsymbol A}}_K^{\rm T}{\boldsymbol P}\bar{{\boldsymbol I}}{\boldsymbol w}(k)+{\boldsymbol w}^{\rm T}(k)\bar{{\boldsymbol I}}^{\rm T}{\boldsymbol P}\bar{{\boldsymbol I}}{\boldsymbol w}(k)+c\big) =\\&\qquad{\cal{E}} \Big( \bar{{\boldsymbol x}}^{\rm T}(k){\left(\bar{{\boldsymbol C}}_K^{\rm T}{\boldsymbol Q}\bar{{\boldsymbol C}}_K+\bar{{\boldsymbol K}}^{\rm T}{\boldsymbol R}\bar{{\boldsymbol K}}+\gamma \bar{{\boldsymbol A}}_K^{\rm T} {\boldsymbol P} \bar{{\boldsymbol A}}_K\right)}\Big.\times\;\\&\qquad\Big.{\bar{{\boldsymbol x}}(k)}\Big)+{\rm{tr}}({\boldsymbol \Sigma}_\eta {\boldsymbol Q}) +\gamma {\rm{tr}}({\boldsymbol \Sigma}_w \bar{{\boldsymbol I}}^{\rm T}{\boldsymbol P}\bar{{\boldsymbol I}})+\gamma c \end{split} $$ □ 由定理1可知, 指标
$$ V(k) = {\cal{E}} {\bar{{\boldsymbol x}}^{\rm T}(k)}{\boldsymbol P}{\bar{{\boldsymbol x}}(k)}+c $$ 可预测当前控制器在时间段$ [k,\infty) $内的控制性能. 当系统的模型可以精确获取时, 可通过求解式(15)和式(16)来计算系统控制性能预测矩阵$ {\boldsymbol P} $和$ c $.
当系统参数发生变化或控制器参数不匹配时, 相应的控制性能和对应的性能预测矩阵$ {\boldsymbol P} $也会随之改变. 由于变化后的系统模型未知, 使得基于模型的计算方法不再适用于$ {\boldsymbol P} $的计算. 随着控制系统和各种智能化仪表及现场总线技术在工业过程中的广泛应用, 大量的过程数据被采集并存储下来, 因此本文在不辨识系统模型的前提下利用采集的过程数据对系统性能预测矩阵进行直接辨识.
为了达到上述目的, 首先令
$$ \bar{{\boldsymbol x}}_s(k)=\left[ {\begin{array}{*{20}{c}} {\boldsymbol x}(k)\\ {\boldsymbol x}_v(k) \end{array}} \right] $$ (21) 则闭环系统的状态空间方程可表示为
$$ \bar{{\boldsymbol x}}_s(k+1) = \bar{{\boldsymbol A}}\bar{{\boldsymbol x}}_s(k)+\bar{{\boldsymbol B}}{\boldsymbol r}(k)+\bar{{\boldsymbol I}}_1{\boldsymbol w}(k) $$ (22) $$ {\boldsymbol y}(k) = \bar{{\boldsymbol C}}\bar{{\boldsymbol x}}_s(k)+{\boldsymbol\eta}(k)\quad\quad\quad\qquad\qquad\; $$ (23) $$ {\boldsymbol u}(k) = \bar{{\boldsymbol K}}_1\bar{{\boldsymbol x}}_s(k) \qquad\qquad\qquad\qquad\qquad $$ (24) 其中
$$\begin{split} &\bar{{\boldsymbol A}}=\left[ {\begin{array}{*{20}{c}} {\boldsymbol A}+{\boldsymbol B}{\boldsymbol K} & {\boldsymbol B}{\boldsymbol C}_v \\ {\bf 0} & {\boldsymbol A}_v \end{array}} \right],\bar{{\boldsymbol B}}=\left[ {\begin{array}{*{20}{c}} {\bf 0} \\ {\boldsymbol B}_v \end{array}} \right],\bar{{\boldsymbol I}}_1=\left[ {\begin{array}{*{20}{c}} {\boldsymbol I} \\ {\bf 0} \end{array}} \right]\notag\\ &\bar{{\boldsymbol C}}=\left[ {\begin{array}{*{20}{c}} {\boldsymbol C}+{\boldsymbol D}{\boldsymbol K} & {\boldsymbol D}{\boldsymbol C}_v \end{array}} \right],\bar{{\boldsymbol K}}_1=\left[ {\begin{array}{*{20}{c}} {\boldsymbol K}&{\boldsymbol C}_v \end{array}} \right]\notag \end{split}$$ 考虑到性能预测矩阵$ {\boldsymbol P} $可以等效地表示为
$$\begin{array}{*{20}{c}} {\boldsymbol P}=\left[ {\begin{array}{*{20}{c}} {\boldsymbol P}_0 & {\boldsymbol P}_1\\{\boldsymbol P}_1^{\rm T} & {\boldsymbol P}_2 \end{array}} \right] \end{array}$$ 其中, $ {\boldsymbol P}_0\in {\bf R}^{(n+k_v)\times (n+k_v)} $, $ {\boldsymbol P}_1\in {\bf R}^{(n+k_v)\times k_r} $, ${\boldsymbol P}_2\in {\bf R}^{k_r\times k_r}$, 则性能预测指标(14)等价于
$$ \begin{split} V(k) =\;& {\cal{E}} \left({\bar{{\boldsymbol x}}_s^{{\rm{T}}}}(k){\boldsymbol P}_{0}\bar{{\boldsymbol x}}_s(k)+2{\bar{{\boldsymbol x}}_s^{{\rm{T}}}}(k){\boldsymbol P}_{1}{\boldsymbol r}(k)\right)+\\ &{\boldsymbol r}^{\rm T}(k){\boldsymbol P}_{2}{\boldsymbol r}(k)+c \end{split} $$ 由此可知
$$ \begin{split} &{\cal{E}} \left({\boldsymbol y}^{\rm T}(k){\boldsymbol Q}{\boldsymbol y}(k)+{\boldsymbol u}^{\rm T}(k){\boldsymbol R}{\boldsymbol u}(k)\right)+\gamma V(k+1)= \\&\quad\;{\cal{E}}\left( \bar{{\boldsymbol x}}^{\rm T}_s(k)\left(\bar{{\boldsymbol C}}^{\rm T}{\boldsymbol Q}\bar{{\boldsymbol C}}+\bar{{\boldsymbol K}}_1^{\rm T}{\boldsymbol R}\bar{{\boldsymbol K}}_1\right){\bar{{\boldsymbol x}}_s(k)}\;+\right.\\&\quad\;\left.{\boldsymbol\eta}^{\rm T}(k){\boldsymbol Q}{\boldsymbol\eta}(k)\right)+\\&\quad\;\gamma {\cal{E}} \left( \bar{{\boldsymbol x}}^{\rm T}_s(k)\bar{{\boldsymbol A}}^{\rm T} {\boldsymbol P}_{0} \bar{{\boldsymbol A}} {\bar{{\boldsymbol x}}_s(k)}+2\bar{{\boldsymbol x}}^{\rm T}_s(k)\bar{{\boldsymbol A}}^{\rm T}{\boldsymbol P}_{0}\bar{{\boldsymbol B}}{\boldsymbol r}(k)\;+\right. \\&\quad\;{\boldsymbol w}^{\rm T}(k)\bar{{\boldsymbol I}}_1^{\rm T}{\boldsymbol P}_{0}\bar{{\boldsymbol I}}_1{\boldsymbol w}(k)+{\boldsymbol r}^{\rm T}(k)\bar{{\boldsymbol B}}^{\rm T}{\boldsymbol P}_{0}\bar{{\boldsymbol B}}{\boldsymbol r}(k)\;+\\&\quad\;2\bar{{\boldsymbol x}}_s^{\rm T}(k)\bar{{\boldsymbol A}}^{\rm T}{\boldsymbol P}_{1}{\boldsymbol r}(k)+2{\boldsymbol r}^{\rm T}(k)\bar{{\boldsymbol B}}^{\rm T}{\boldsymbol P}_{1}{\boldsymbol r}(k)\;+\\&\quad\;{\boldsymbol r}^{\rm T}(k){\boldsymbol P}_{2}{\boldsymbol r}(k)+2\bar{{\boldsymbol x}}^{\rm T}_s(k){\boldsymbol P}_{0}\bar{{\boldsymbol I}}_1{\boldsymbol w}(k)\;+\\&\quad\left.2{\boldsymbol r}^{\rm T}(k)\bar{{\boldsymbol B}}^{\rm T}{\boldsymbol P}_{0}\bar{{\boldsymbol I}}_1{\boldsymbol w}(k)+c\right) \end{split} $$ 综合Bellman方程, 可得
$$ {\boldsymbol P}_{0} = \bar{{\boldsymbol C}}^{\rm T}{\boldsymbol Q}\bar{{\boldsymbol C}}+\bar{{\boldsymbol K}}_1^{\rm T}{\boldsymbol R}\bar{{\boldsymbol K}}_1+\gamma \bar{{\boldsymbol A}}^{\rm T}{\boldsymbol P}_{0}\bar{{\boldsymbol A}} $$ (25) $$ {\boldsymbol P}_{1} = \gamma \bar{{\boldsymbol A}}^{\rm T}{\boldsymbol P}_{0}\bar{{\boldsymbol B}}+\gamma \bar{{\boldsymbol A}}^{\rm T}{\boldsymbol P}_{1}\qquad\qquad\quad\; $$ (26) $$ {\boldsymbol P}_{2} = \gamma {\boldsymbol P}_{2}+2\gamma \bar{{\boldsymbol B}}^{\rm T}{\boldsymbol P}_{1} +\gamma \bar{{\boldsymbol B}}^{\rm T}{\boldsymbol P}_{0}\bar{{\boldsymbol B}}\quad\quad $$ (27) 由式(27)可得
$$ {\boldsymbol P}_{2} = \frac{2\gamma \bar{{\boldsymbol B}}^{\rm T}{\boldsymbol P}_{1} +\gamma \bar{{\boldsymbol B}}^{\rm T}{\boldsymbol P}_{0}\bar{{\boldsymbol B}}}{1-\gamma} $$ (28) 显然, $ {\boldsymbol P}_{2} $可由参数$ {\boldsymbol P}_{0},{\boldsymbol P}_{1} $求得. 因此, 只需要辨识参数$ {\boldsymbol P}_{0},{\boldsymbol P}_{1} $即可得到$ {\boldsymbol P}_2 $和性能预测矩阵$ {\boldsymbol P} $, 这样可以减少要辨识参数的个数, 从而提高辨识精度.
为了采用数据驱动方法辨识$ {\boldsymbol P}_{0},{\boldsymbol P}_{1} $, 用$ \bar{{\boldsymbol x}}_s(k) $代替其均值, 则
$$ \begin{split} V&(k) = {\bar{{\boldsymbol x}}_s^{{\rm{T}}}}(k){\boldsymbol P}_{0}\bar{{\boldsymbol x}}_s(k)+2{\bar{{\boldsymbol x}}_s^{{\rm{T}}}}(k){\boldsymbol P}_{1}{\boldsymbol r}(k)\notag\;+\\ &{\boldsymbol r}^{\rm T}(k)\frac{2\gamma \bar{{\boldsymbol B}}^{\rm T}{\boldsymbol P}_{1} +\gamma \bar{{\boldsymbol B}}^{\rm T}{\boldsymbol P}_{0}\bar{{\boldsymbol B}}}{1-\gamma}{\boldsymbol r}(k)+c \notag=\\ &{\bar{{\boldsymbol x}}_s^{{\rm{T}}}}(k){\boldsymbol P}_{0}\bar{{\boldsymbol x}}_s(k)+2\left({\bar{{\boldsymbol x}}_s}(k)+\frac{\gamma}{1-\gamma} \bar{{\boldsymbol B}} {\boldsymbol r}(k)\right)^{\rm T}\times\\ &{\boldsymbol P}_{1}{\boldsymbol r}(k)\notag\;+{\boldsymbol r}^{\rm T}(k)\frac{\gamma \bar{{\boldsymbol B}}^{\rm T}{\boldsymbol P}_{0}\bar{{\boldsymbol B}}}{1-\gamma}{\boldsymbol r}(k)+c \notag=\\ &(\bar{{\boldsymbol x}}_s^{\rm T}(k)\otimes \bar{{\boldsymbol x}}_s^{\rm T}(k)){\boldsymbol D}_{n+k_v}hvec({\boldsymbol P}_{0})\notag\;+\\ &2\left(\left({\bar{{\boldsymbol x}}_s}(k)+\frac{\gamma}{1-\gamma} \bar{{\boldsymbol B}} {\boldsymbol r}(k)\right)^{\rm T} \otimes {\boldsymbol r}^{\rm T}(k)\right)vec({\boldsymbol P}_{1})\notag\;+\\ &\frac{\gamma}{1-\gamma}\left({\boldsymbol r}^{\rm T}(k)\bar{{\boldsymbol B}}^{\rm T}\otimes{\boldsymbol r}^{\rm T}(k)\bar{{\boldsymbol B}}^{\rm T}\right){\boldsymbol D}_{n+k_v}hvec({\boldsymbol P}_{0})+c \notag \end{split}$$ 其中, $ \otimes $表示Kronecker积. $ hvec({\boldsymbol P}_{0}) $表示由对称矩阵$ {\boldsymbol P}_{0} $的上三角矩阵的列叠加而成的向量, 表示待辨识的$(n+k_v)(n+k_v+1)/{2}$个参数, 且满足
$$ \begin{split} &{\boldsymbol D}_{n+k_v}hvec({\boldsymbol P}_0) = vec({\boldsymbol P}_0)\\ &{\boldsymbol D}_{n+k_v}\in {\bf R}^{(n+k_v)^2\times\tfrac{(n+k_v)(n+k_v+1)}{2}} \end{split} $$ $ {\boldsymbol D}_{n+k_v} $称为重复矩阵[39]. $ vec({\boldsymbol P}_{1}) $表示由对称矩阵$ {\boldsymbol P}_{1} $的列叠加而成的向量. 令
$$ \begin{split} &{\boldsymbol p}_1 = hvec({\boldsymbol P}_{0}), \;\;{\boldsymbol p}_2 = vec({\boldsymbol P}_{1}) \\ &{\boldsymbol \phi}_1(k) = \Bigg(\bar{{\boldsymbol x}}_s^{\rm T}(k)\otimes\bar{{\boldsymbol x}}^{\rm T}_s(k)\;+\\ &\;\;\quad\quad\quad\frac{\gamma}{1-\gamma}{\boldsymbol r}^{\rm T}(k) \bar{{\boldsymbol B}}^{\rm T}\otimes {\boldsymbol r}^{\rm T}(k)\bar{{\boldsymbol B}}^{\rm T}\Bigg){\boldsymbol D}_{n+k_v} \\ &{\boldsymbol \phi}_2(k) = 2\left({\bar{{\boldsymbol x}}_s}(k)+\frac{\gamma}{1-\gamma} \bar{{\boldsymbol B}} {\boldsymbol r}(k)\right)^{\rm T}\otimes {\boldsymbol r}^{\rm T}(k) \end{split} $$ 则Bellman方程可等效表示为
$$ \begin{split} {{\boldsymbol \phi}_1(k)}&{{\boldsymbol p}_1} + {\boldsymbol \phi}_2(k){\boldsymbol p}_2+c = \\ &{\boldsymbol u}^{\rm T}(k){\boldsymbol R}{\boldsymbol u}(k)+{\boldsymbol y}^{\rm T}(k){\boldsymbol Q}{\boldsymbol y}(k)\;+\\ &\gamma ({\boldsymbol \phi}_1(k+1){\boldsymbol p}_1+{\boldsymbol \phi}_2(k+1){\boldsymbol p}_2+c) \end{split} $$ 上式等价于
$$ \begin{split} ({\boldsymbol \phi}_1(k)- & \gamma{\boldsymbol \phi}_1(k+1)){\boldsymbol p}_1 +({\boldsymbol \phi}_2(k)- \gamma{\boldsymbol \phi}_2(k+1)){\boldsymbol p}_2\;+\\ &(1-\gamma)c = {\boldsymbol u}^{\rm T}(k){\boldsymbol R}{\boldsymbol u}(k)+{\boldsymbol y}^{\rm T}(k){\boldsymbol Q}{\boldsymbol y}(k)) \end{split} $$ 令
$$ \begin{split} &{\boldsymbol \theta} = \left[ {\begin{array}{*{20}{c}} {\boldsymbol{p}}_1^{\rm{T}} & {\boldsymbol{p}}_2^{\rm{T}} &c \end{array}} \right]^{\rm T}\\ &{\boldsymbol {\varphi}}(k) =\\ & \qquad\ \left[ {\boldsymbol{\phi}}_1(k)-\lambda{\boldsymbol{\phi}}_1(k+1)\;\; {\boldsymbol{\phi}}_2(k)-\lambda{\boldsymbol{\phi}}_2(k+1) \;\;1-\lambda \right]\\ &{\boldsymbol \kappa}(k) = {\boldsymbol u}^{\rm T}(k){\boldsymbol R}{\boldsymbol u}(k)+{\boldsymbol y}^{\rm T}(k){\boldsymbol Q}{\boldsymbol y}(k) \end{split} $$ 可得
$$ {\boldsymbol \varphi}(k){\boldsymbol \theta} = {\boldsymbol \kappa}(k) $$ 则$ {\boldsymbol P} $可由算法1求得. 实际应用中也可借助递推最小二乘法利用在线采集的数据对$ {\boldsymbol \theta} $和矩阵$ {\boldsymbol P} $进行实时更新. 总的来说, 对反馈控制系统而言, 系统的性能指标可用关键指标$ {\boldsymbol P} $来预测.
算法1. ${\boldsymbol P}$的在线辨识
1) 计算并收集过程数据${\boldsymbol u}(i),{\boldsymbol y}(i),{\boldsymbol x}(i),{\boldsymbol x}_v(i),{\boldsymbol r}(i),i= k_0,k_0+1,\cdots ,\;k_0+N_0$;
2) 形成
$$ \begin{array}{l} {\boldsymbol \Phi}=\left[ {\begin{array}{*{20}{c}} {\boldsymbol \varphi}(k_0)\\ {\boldsymbol \varphi}(k_0+1)\\ \vdots\\ {\boldsymbol \varphi}(k_0+N_0) \end{array}} \right],\;\;{\boldsymbol \Gamma}=\left[ {\begin{array}{*{20}{c}} {\boldsymbol \kappa}(k_0)\\ {\boldsymbol \kappa}(k_0+1)\\ \vdots\\ {\boldsymbol \kappa}(k_0+N_0) \end{array}} \right] \end{array} $$ 3) 求解
$${\boldsymbol \Phi}{\boldsymbol \theta} = {\boldsymbol \Gamma} \Longrightarrow \hat{{\boldsymbol \theta}} = \left({\boldsymbol \Phi}^{\rm T}{\boldsymbol \Phi} \right)^{-1}{\boldsymbol \Phi}^{\rm T} {\boldsymbol \Gamma} $$ 4) 根据$ \hat{{\boldsymbol \theta}} $构建矩阵$ {\boldsymbol P}_{0},{\boldsymbol P}_{1} $;
5) 利用式(28)计算$ {\boldsymbol P}_{2} $;
6) 构建性能预测矩阵${\boldsymbol P}=\left[ \begin{aligned}{\boldsymbol P}_0 \;\;\;{\boldsymbol P}_1\\{\boldsymbol P}_1^{\rm T} \;\;{\boldsymbol P}_2 \end{aligned}\right]$.
2.2 基于黎曼度量的控制性能监测
值得注意的是, 系统参数矩阵$ {\boldsymbol A}, {\boldsymbol B}, {\boldsymbol C}, {\boldsymbol D} $与控制器参数的改变都会引起性能预测矩阵$ {\boldsymbol P}_0 $的变化. 为了更好地阐述$ {\boldsymbol P}_0 $的变化与系统动态特性变化之间的关系, 首先考虑Lyapunov等式(25)的一种简单形式
$$\left\{ \begin{aligned} &\bar{\boldsymbol A}^{\rm T}{\boldsymbol P}_0\bar{\boldsymbol A}-{\boldsymbol P}_0+\bar{\boldsymbol W} = 0\\ & \bar{\boldsymbol W} = \bar{\boldsymbol C}^{\rm T}{\boldsymbol Q}\bar{\boldsymbol C}+\bar{\boldsymbol {{K}}}_1^{\rm T}{\boldsymbol R}\bar{{\boldsymbol K}}_1 \end{aligned}\right. $$ (29) 式(29)表示了系统无噪声时$ \gamma $选为$ 1 $的特殊情况. 由文献[38]可知, 矩阵$ {\boldsymbol P}_0 $可以表示为如下形式
$$ {\boldsymbol P}_0 = \sum\limits_{i = k}^\infty (\bar{{\boldsymbol A}}^i)^{\rm T}\bar{{\boldsymbol W}}\bar{{\boldsymbol A}}^i $$ (30) 显然, 系统矩阵$ {\boldsymbol A},{\boldsymbol B},{\boldsymbol C},{\boldsymbol D} $与控制器参数的变化都体现在性能预测矩阵$ {\boldsymbol P}_0 $的变化中. 因此, 本文主要关注矩阵$ {\boldsymbol P}_0 $的变化, 并借助矩阵$ {\boldsymbol P}_0 $来实现系统控制性能监测与诊断.
由上述可知, 性能预测矩阵$ {\boldsymbol P}_0 $包含了系统在正常操作状态、故障模态、不确定性等情况下的系统性能衰退与改变的信息. 由于该矩阵$ {\boldsymbol P}_0\in {\bf R}^{n\times n} $具有对称、正定的结构属性, 这种特殊结构属于正定对称矩阵的${n\times(n+1)}/{2}$维黎曼流形. 由于黎曼度量可以有效地衡量流形曲面上两个正定矩阵在方向和幅值方面的变化, 因此本文利用黎曼度量来衡量$ {\boldsymbol P}_0 $矩阵的实时变化, 进而实现系统控制性能监测与诊断.
由于系统存在不确定性, 不同的不确定性参数对应不同的性能预测矩阵. 为了利用黎曼度量来实现性能监测, 首先需要确认含有随机不确定性的系统中心(黎曼均值)并设定阈值. 从式(9)和式(10)可知, 这需要知道所有可能的$ {\boldsymbol P}_0 $才能实现. 然而, 由于参数$ \delta_i $在区间 $ [\bar{\delta}_{i,1},\bar{\delta}_{i,2}] $内变化, 使得$ {\boldsymbol P}_0 $有无数种可能, 为阈值的设定带来困难. 近年来, 随机算法(Randomized algorithm, RA)为随机框架下极值的估计提供了有效的工具[40-41]. 因此, 本文基于RA算法来实现阈值设定. 为了上述目的, 首先引入定理 2.
定理 2 [40, 42]. 给定$ \alpha\in (0,1) $, $ \epsilon\in (0,1) $, 令
$$ N\ge \frac{\ln \frac{1}{\epsilon}}{\ln \frac{1}{1-\alpha}} $$ (31) 根据随机变量$ w $的概率密度函数$ D(w) $和支持度$ {\cal{D}}_w $生成$ N $个独立同分布的变量$ w_i, i = 1,\cdots ,\;N $, 则
$$ \begin{split} &{\rm{Pr}}(\Theta(w)\le \gamma_{\max})\ge 1-\alpha \\ &\gamma_{\max} = \max_{i = 1,\cdots ,\;N} \Theta(w_i)\notag \end{split} $$ 的置信度大于$ 1-\epsilon $, 其中, $ \Theta(w) $是关于随机变量$ w $的函数.
令
$$ \bar{\delta}=\left[ {\begin{array}{*{20}{c}} \delta_1\\ \vdots\\ \delta_l \end{array}} \right],\;\;\bar{\delta}_1=\left[ {\begin{array}{*{20}{c}} \bar{\delta}_{1,1}\\ \vdots\\ \bar{\delta}_{l,1} \end{array}} \right],\;\;\bar{\delta}_2=\left[ {\begin{array}{*{20}{c}} \bar{\delta}_{1,2}\\ \vdots\\ \bar{\delta}_{l,2} \end{array}} \right] $$ (32) 则$ \delta $为区间 $ [\bar{\delta}_1,\bar{\delta}_2] $内分布的随机变量. 由定理2可知, 算法2可用于实现阈值$ J_{{\rm{th}}} $的设定.
算法2. 基于RA算法的阈值设定
1) 根据给定$ \alpha\in (0,1) $, $ \epsilon\in (0,1) $, 利用
$$ {N\ge\frac{\ln \frac{1}{\epsilon}}{\ln \frac{1}{1-\alpha}}} $$ 确定整数$ N $;
2) 生成$ N $个随机样本$\bar{\delta}_i\sim {\rm{U}}[\bar{\delta}_1,\bar{\delta}_2],i=1,\cdots ,\; N$, 其中$ {\rm{U}}[\bar{\delta}_1,\bar{\delta}_2] $表示区间$ [\bar{\delta}_1,\bar{\delta}_2] $内的均匀分布;
3) 根据$ \bar{\delta}_i $生成对应的$ {\boldsymbol A}(\bar{\delta}_i), {\boldsymbol B}(\bar{\delta}_i), {\boldsymbol C}(\bar{\delta}_i), {\boldsymbol D}(\bar{\delta}_i) $, 根据式(15)计算并提取对应的$ {\boldsymbol P}_0(\bar{\delta}_i), i=1,\cdots ,\;N $;
4) 确定黎曼均值
$$ {\boldsymbol P}_z=\arg \min\limits_{{\boldsymbol P}_j, j=1,\cdots ,\; N} \sum\limits_{i=1}^N d_R^2({\boldsymbol P}_i, {\boldsymbol P}_j) $$ 5) 设定阈值
$$ J_{{\rm{th}}}=\max\limits_{i\in \{1,\cdots ,\;N\}} d_R^2({\boldsymbol P}_i,{\boldsymbol P}_z) $$ 算法2保证在无故障情况下$, $ 评估函数
$$ J = d_R^2({\boldsymbol P}_0,{\boldsymbol P}_z) $$ (33) 大于阈值的概率
$$ {\rm{Pr}}(J>J_{{\rm{th}}}|{\text{无故障}})\ge 1-\alpha $$ (34) 的置信度大于$ 1-\epsilon $.
在故障诊断中,
$$ p_{{\rm{FAR}}} = {\rm{Pr}}(J>J_{{\rm{th}}}|{\text{无故障}}) $$ (35) 通常称为误报率(False alarm rate, FAR)[2]. 也就是说, 通过算法2实现的性能监测方法的误报率小于$ \alpha $的置信度大于$ 1-\epsilon $.
在性能预测与RA算法的基础上, 本文提出了基于黎曼度量的控制性能监测方案. 该方案分为离线建模与在线检测两部分.
离线建模过程为:
1)利用RA算法生成$ N $个随机不确定性系统模型(1);
2)利用算法2确定黎曼均值$ {\boldsymbol P}_z $和阈值${J_{{\rm{th}}}}$.
在线检测过程为:
1)采集现场数据, 通过算法1对实时的系统性能预测矩阵${\boldsymbol P}_{{\rm{new}}}$进行辨识;
2)计算实时性能矩阵${{\boldsymbol P}_{{\rm{new}}}}$与黎曼均值$ {{\boldsymbol P}_z} $的黎曼度量, 并将其作为检测指标$ J $
$$ \begin{array}{l} J = d_R^2({\boldsymbol P}_{{\rm{new}}},{\boldsymbol P}_z) \end{array} $$ (36) 3)运行性能监测逻辑
$$ \begin{array}{l} \begin{cases} J_{{\rm{th}}}-J> 0 \Longrightarrow {\text{系统无异常}}\\ J_{{\rm{th}}}-J\le 0 \Longrightarrow {\text{系统性能衰退}}\end{cases} \end{array} $$ (37) 总结而言, 基于黎曼度量的控制性能监测流程如图2所示.
图 2 基于黎曼度量的控制性能监测流程图Fig. 2 Flow chart of Riemannian metric based control performance monitoring值得注意的是, 本文所提出的基于黎曼度量的控制性能监测方法无需辨识发生性能衰退后系统的模型, 仅通过在线识别其性能预测矩阵${\boldsymbol P}_{{\rm{new}}}$就能实现对系统性能衰退程度的监测. 尽管系统性能预测矩阵${\boldsymbol P}_{{\rm{new}}}$的辨识需要用到状态变量, 但${\boldsymbol P}_{{\rm{new}}}$仅取决于系统与控制器参数的变化, 对状态变量的变化具有不变性. 另外, 由于Bellman方程的在线求解对控制策略的在线优化至关重要[43], 因此该性能预测矩阵也可以用于有效地指导容错控制器的设计和优化.
在文献[42]中, Ding等利用随机算法给出了故障诊断系统的阈值设定与诊断性能评估方法, 该方法主要借助基于观测器的残差发生器所生成的残差信号进行故障诊断. 而本文主要针对带有反馈控制环节的动态系统, 提出一种系统性能的实时预测指标, 并通过黎曼度量衡量实时性能预测指标的退化程度来实现系统性能监测与诊断. 与文献[42]相比, 本文具有如下的创新性:
1)提出了反馈控制系统的性能衰退预测方法并给出了对应的在线辨识方案;
2)利用黎曼度量实现了控制性能监测与诊断.
3. 基于黎曼度量的故障隔离
当检测到系统发生性能衰退后, 紧接着需要进行的就是对引发系统性能衰退的故障进行隔离, 判断当前性能衰退的类型, 并根据故障类型及时采取相应措施. 对于加性故障而言, 通常根据不同的故障类型设计对应的残差发生器, 然后利用残差评估函数和决策逻辑实现故障分离. 这类方法的应用关键在于: 1)对要分离的故障进行聚类; 2)设计一系列残差发生器保证每一个残差发生器只对其中一类故障敏感; 3)针对每一个残差发生器设计对应的阈值[5]. 目前对乘性故障隔离的研究相对较少, 因此本节设计基于黎曼度量的乘性故障的隔离方法.
假设系统有$ M $个故障模态, 对于每一个类型的故障对应一个性能预测矩阵的集合$\mathcal P_i,i = 1,\cdots ,\; M$. 该集合中的性能预测矩阵可通过故障和不确定性模型求得. 为了实现故障分离, 首先需要求取每一个集合$ \mathcal P_i $的黎曼均值和对应的阈值. 假设集合$ \mathcal P_i $中性能预测矩阵为$ {\boldsymbol P}_{i,n},n = 1,\cdots ,\;N_i $, 其中$ N_i $为集合$ {\cal{P}}_i $中的性能预测矩阵个数. 则其黎曼均值可通过求解式(38)得到
$$ {\boldsymbol P}_{z,i} = \arg \min_{{\boldsymbol P}_{i,j}\in{\cal{P}}_i} \sum\limits_{n = 1}^{N_i} d_R^2({\boldsymbol P}_{i,j},{\boldsymbol P}_{i,n}) $$ (38) 本文将这个黎曼均值称为该故障模态的中心. 所有矩阵与故障模态中心的黎曼度量中的最大值为
$$ \gamma_{i} = \max\limits_{{\boldsymbol P}_{i,j}\in{\cal{P}}_i} d_R^2({\boldsymbol P}_{i,j},{\boldsymbol P}_{z,i}) $$ (39) 则以$ {\boldsymbol P}_{z,i} $为中心、$ \gamma_i $为半径, 可定义如下故障簇
$$ {\cal{C}}_{f_i} = \left\{ {\boldsymbol P}:d_R^2({\boldsymbol P},{\boldsymbol P}_{z,i})\le \gamma_i\right\} $$ (40) 本文中假设任一故障仅属于一个故障簇. 当在线检测到故障时, 可通过辨识的故障系统的性能预测矩阵$ {\boldsymbol P}_{{\rm{new}}} $, 利用如下决策逻辑实现故障隔离
$$ \begin{array}{l} \begin{cases} d_R^2({\boldsymbol P}_{{\rm{new}}},{\boldsymbol P}_{z,i})\le \gamma_{i}\\ d_R^2({\boldsymbol P}_{{\rm{new}}},{\boldsymbol P}_{z,j})> \gamma_{j},\;\;j\neq i,j\in [1,M] \end{cases} \end{array} $$ 等价于故障$ i $发生.
4. 三容水箱系统的性能监测与诊断
本节利用三容水箱的实验平台(DTS200)对所提出的算法进行验证.
4.1 实验平台介绍
三容水箱具有化工过程中常用到的储罐、管道和泵, 是典型的过程控制实验设备, 如图3所示. 三容水箱既可以用来模拟工业生产过程中的液位控制, 也能够模拟各种实际应用中的典型故障, 如传感器失效、执行器失效、水箱漏水、连通阀阻塞等, 因此在故障检测研究中也得到了广泛应用.
三容水箱可用如下数学模型描述
$$ \begin{split} &\dot{x}_{1} = {\cal{A}}_{\mathrm{inv}}u_1-{\cal{A}}_{\mathrm{inv}}a_{1}s_{13}\text{sgn}(x_{1}-x_{3})\sqrt{2g|x_{1}-x_{3}|}\\ &\dot{x}_{2} = {\cal{A}}_{\mathrm{inv}}u_{2}+{\cal{A}}_{\mathrm{inv}}a_{3}s_{23}\text{sgn}(x_{3}-x_{2}) \sqrt{2g|x_{3}-x_{2}|}\;-\\ &\;\qquad{\cal{A}}_{\mathrm{inv}}s_{0}\sqrt{2gx_{2}}\\ &\dot{x}_{3} = {\cal{A}}_{\mathrm{inv}}a_{1}s_{13}\text{sgn}(x_{1}-x_{3}) \sqrt{2g|x_{1}-x_{3}|}\;-\\ &\;\qquad{\cal{A}}_{\mathrm{inv}}a_{3}s_{23}\text{sgn}(x_{3}-x_{2})\sqrt{2g|x_{3}-x_{2}|} \end{split} $$ 其中, $ x_{i}(t) = h_i(t),i = 1,2,3 $表示每个水箱的水位, ${\cal{A}} $表示水箱面积, $ {\cal{A}}_{\mathrm{inv}} = 1/{\cal{A}} $且$ s_{13} = s_{23} = s_{0} = s_{n} $. $ u_1,u_2 $分别表示水泵1和水泵2的进水量. 模型参数通过实际实验平台测量得到, 如表1所示. 三容水箱在操作点$h_1 = 25\;{\rm{cm}}$, $h_2 = 20\;{\rm{cm}}$, $h_3 = 22.5\;{\rm{cm}}$附近工作, 水箱数据的采样周期为$2\,{\rm{s}}$. 此时水箱的系统模型可用式(1)进行描述, 其中
表 1 水箱DTS200的参数Table 1 Parameters of tank DTS200参数 符号 值 单位 水箱面积 $ {\cal{A}} $ 154 $ \mathrm{cm}^{2} $ 水箱间管道面积 $s_{{n} }$ 0.5 $ \mathrm{cm}^{2} $ 水箱最高水位 $ H_{\mathrm{max}} $ 62 $ \mathrm{cm} $ 泵 1 的最大进水量 $ Q_{1_{\mathrm{max}}} $ 120 $ \mathrm{cm}^{3}/\mathrm{s} $ 泵 2 的最大进水量 $ Q_{2_{\mathrm{max}}} $ 120 $ \mathrm{cm}^{3}/\mathrm{s} $ 管道 1 水流系数 $ a_{1} $ 0.45 管道 2 水流系数 $ a_{2} $ 0.60 管道 3 水流系数 $ a_{3} $ 0.45 $$\begin{split} & {\boldsymbol A}_0=\left[ {\begin{array}{*{20}{c}} 0.9591 & -0.0005 & 0.0370\\ -0.0037 & 0.9388 & 0.0341\\ 0.0384 & 0.0380 & 0.9229 \end{array}} \right],\;{\boldsymbol D}_0=0 \\ & {\boldsymbol B}_0=\left[ {\begin{array}{*{20}{c}} 0.0127 & 0\\ 0 & 0.0127\\ 0 & 0 \end{array}} \right],\; {\boldsymbol C}_0=\left[ {\begin{array}{*{20}{c}} 1 & 0 & 0\\0 & 1 & 0\\0 & 0 &1 \end{array}} \right] \end{split} $$ 这里, 采用如下反馈控制器
$$ {\boldsymbol K}=-\left[ {\begin{array}{*{20}{c}} 0.9329 & 0.0048 &0.0982\\ 0.0049 & 0.8876 & 0.0967 \end{array}} \right] $$ 由于水箱参数$ a_1,a_3 $往往无法精确确定, 可能在一定范围内波动, 如$ a_1\in[0.35,0.55] $, $a_3\in [0.35, 0.55]$. 因此不可避免地会给线性模型带来不确定性, 这个不确定性可以用凸多面体不确定性来描述, 其中
$$\begin{split} &{{\boldsymbol \Delta} _{{A_1}}}=\left[ {\begin{array}{*{20}{c}} 0.0088 & 0.0001 & -0.0079\\ -0.0002 & 0 & 0.0002\\ -0.0086 & -0.0001 & 0.0079 \end{array}} \right] \\ & {{\boldsymbol \Delta} _{{A_2}}}=\left[ {\begin{array}{*{20}{c}} -0.0087 & -0.0001 & 0.0079\\ 0.0002 & 0 & -0.0002\\ 0.0085 & 0.0001 & 0.0077 \end{array}} \right] \\ &{{\boldsymbol \Delta} _{{A_3}}}=\left[ {\begin{array}{*{20}{c}} 0 & -0.0002 & 0.0002\\ -0.0005 & 0.0080 & -0.0085\\ 0.0005 & -0.0079 & 0.0084 \end{array}} \right] \\ &{{\boldsymbol \Delta} _{{A_4}}}=\left[ {\begin{array}{*{20}{c}} 0 & 0.0002 & -0.0001\\ 0.0006 & -0.0079 & 0.0084\\ -0.0005 & 0.0078 & -0.0084 \end{array}} \right]\end{split}$$ 4.2 性能监测结果
对三容水箱而言, 液位高度是衡量水箱运行状态的关键. 因此在性能指标中, 液位高度的加权比重选得比较大, 具体为$ Q = 500,\; R = 1 $. 给定故障误报率$ \alpha = 0.01 $, 置信度的显著性水平$ \delta = {10^{ - 7}} $, 利用RA算法可得$ N $应该不小于$1\,604$. 运行算法2可得阈值$ {J_{{\rm{th}}} = 0.0524} $, 对应的黎曼均值为
$$ \begin{array}{l} {\boldsymbol P}_z=10^3\times \left[ {\begin{array}{*{20}{c}} 4.4323 & 0.5024 & 1.5905\\ 0.5024 & 3.5168 & 1.2925\\ 1.5905 & 1.2925 & 3.5264 \end{array}} \right] \end{array} $$ 在水箱的长期操作中, 由于存在水垢等的影响, 管道可能会具有一定程度的堵塞. 管道的堵塞是导致三容水箱性能衰退的其中一个关键因素. 因此, 本文主要针对三容水箱的管道堵塞引发的控制性能衰退展开性能监测与诊断的实验验证. 在第$ 2\,000 $个采样点, 管道1发生堵塞, 从而导致水流系数变为$ a_1 = 0.25 $. 采集过程数据, 利用迭代最小二乘算法(算法1的迭代实现)对性能预测矩阵进行在线辨识. 其对应的性能监测效果如图4所示, 显然故障发生后, 评估函数很快高于阈值. 也就是说, 借助黎曼度量可以及时诊断出故障. 类似地, 可得到管道2和管道3发生堵塞(水流系数分别变为$a_2 = 0.3$和$ a_3 = 0.25 $)的诊断结果分别如图5和图6所示. 另一方面, 当控制参数发生不匹配时, 如控制器增益在第$2\;000$个采样点变为原来的两倍, 其对应的性能诊断结果如图7所示. 显然, 所提出的方法不仅能诊断故障引起的性能衰退, 也可诊断控制参数变化引发的性能退化.
4.3 比较实验
基于性能衰退预测的故障诊断是近年来针对反馈控制系统的一个新的研究课题, 由文献[16]首次提出. 该方法的优点是将性能衰退预测的信息用于过程监测, 对由参数变化或乘性故障引起的系统性能衰退的实时检测与控制补偿非常有效. 因此, 本节将本文所提出的方法与文献[16]中提出的基于系统性能退化指标(Indicator for the system performance degradation, ISPD)的性能监测方法进行比较. 在基于ISPD的性能监测算法中, $ \bar{{\boldsymbol x}}^{\rm T}(k){\boldsymbol P}_z\bar{{\boldsymbol x}}(k) $用于对标称系统的性能进行预测, $ \bar{{\boldsymbol x}}^{\rm T}(k){\boldsymbol P}_{{\rm{new}}}\bar{{\boldsymbol x}}(k) $用于对实时系统性能进行预测, 而两者的比值
$$ J_{{\rm{ISPD}}} = \frac{\bar{{\boldsymbol x}}^{\rm T}(k){\boldsymbol P}_{{\rm{new}}}\bar{{\boldsymbol x}}(k)}{\bar{{\boldsymbol x}}^{\rm T}(k){\boldsymbol P}_{z}\bar{{\boldsymbol x}}(k)} $$ (41) 作为ISPD用于对系统性能衰退程度进行监测. 在式(41)中, $ {\boldsymbol P}_{{\rm{new}}} $表示利用算法1的迭代算法辨识的实时性能预测矩阵. 通过RA算法设定对应的阈值$ J_{{\rm{ISPD}},{\rm{th}}} $, 则可通过如下检测逻辑实现性能监测.
$$ \begin{cases} J_{{\rm{ISPD}},{\rm{th}}}-J_{{\rm{ISPD}}}>0 \Longrightarrow {\text{系统无异常}}\\ J_{{\rm{ISPD}},{\rm{th}}}-J_{{\rm{ISPD}}}\le 0 \Longrightarrow {\text{系统性能衰退}} \end{cases} $$ 当管道2发生堵塞时(水流系数变为$ a_2 = 0.3 $), 基于ISPD的性能监测结果如图8所示. 与图5比较可知, 基于黎曼度量的性能监测方法对系统性能衰退的检测能力更强.
实际上, 基于ISPD的性能监测方法和本文所提出的方法都基于性能预测指标来实现性能监测. 两者的不同之处在于基于ISPD的方法监测的是二次型性能指标$ \bar{{\boldsymbol x}}^{\rm T}(k){\boldsymbol P}_z\bar{{\boldsymbol x}}(k) $与$\bar{{\boldsymbol x}}^{\rm T}(k){\boldsymbol P}_{{\rm{new}}}\bar{{\boldsymbol x}}(k)$这两个标量的比值, 而本文通过衡量两个性能预测指标矩阵之间的黎曼距离来实现性能监测. 与标量相比, 黎曼距离能更好地反映性能预测矩阵在方向和幅值方面的改变, 因此可以更有效地检测系统性能的变化, 进而实现性能监测与故障隔离. 这也是基于黎曼度量的性能监测方法与基于ISPD方法相比, 具有更高的可检测性的原因.
4.4 故障隔离结果
当检测到系统性能衰退后, 需要对系统性能衰退的类型进行隔离, 判断故障种类, 以便采取相应措施消除故障. 本节将故障分为3 类: 管道1堵塞、管道2堵塞、管道3堵塞. 利用提出的故障隔离算法, 借助RA算法可得到每一个故障簇的中心
$$ \begin{split} &{\boldsymbol P}_{z,1}=10^3\times \left[ {\begin{array}{*{20}{c}} 4.7222 & 0.4175 & 1.4102\\ 0.4175 & 3.5403 & 1.3440\\ 1.4102 & 1.3440 & 3.6406 \end{array}} \right] \\ &{\boldsymbol P}_{z,2}=10^3\times \left[ {\begin{array}{*{20}{c}} 4.4682 & 0.5971 & 1.6472\\ 0.5971 & 3.8132 & 1.4452\\ 1.6472 & 1.4452 & 3.6122 \end{array}} \right] \\ &{\boldsymbol P}_{z,3}=10^3\times \left[ {\begin{array}{*{20}{c}} 4.5140 & 0.4079 & 1.7140\\ 0.4079 & 3.6688 & 1.1249\\ 1.7140 & 1.1249 & 3.7223 \end{array}} \right] \end{split}$$ 对应的故障半径分别为
$$ \begin{array}{l} \gamma_1 = 0.0808,\; \gamma_2 = 0.0414,\; \gamma_3 = 0.0712 \end{array} $$ 本节分别利用5 个堵塞故障$a_1 = 0.30$, $ a_2 = 0.35 $, $ a_3 = 0.28 $, $ a_1 = 0.27 $, $ a_2 = 0.40 $来验证所提出的故障隔离算法. 对应的故障隔离结果如表2所示.
表 2 水箱堵塞故障隔离Table 2 Isolation of pipe plugging故障 $ d_R^2({\boldsymbol P},{\boldsymbol P}_{z,1}) $ $ d_R^2({\boldsymbol P},{\boldsymbol P}_{z,2}) $ $ d_R^2({\boldsymbol P},{\boldsymbol P}_{z,3}) $ 故障隔离 $a_1 = 0.30$ $ 0.0388 $ $ 0.1088 $ $ 0.1193 $ 故障簇 1 $ a_2 = 0.35 $ $ 0.0894 $ $ 0.0270 $ $ 0.0810 $ 故障簇 2 $ a_3 = 0.28 $ $ 0.1258 $ $ 0.1099 $ $ 0.0478 $ 故障簇 3 $ a_1 = 0.27 $ $ 0.0636 $ $ 0.1325 $ $ 0.1409 $ 故障簇 1 $a_2 = 0.40$ $ 0.0809 $ $ 0.0139 $ $ 0.0732 $ 故障簇 2 由表2可以看出, 当管道1发生堵塞时, 可得
$$ \begin{split} &d_R^2({\boldsymbol P},{\boldsymbol P}_{z,1})<\gamma_1\\ & d_R^2({\boldsymbol P},{\boldsymbol P}_{z,2})>\gamma_2\\ &d_R^2({\boldsymbol P},{\boldsymbol P}_{z,3})>\gamma_3 \end{split} $$ 显然, 该故障属于第1 个故障簇. 类似地, 可以对其余4 个故障进行隔离, 从表2可知, 对这4 个故障进行隔离的结果与真实的故障类型相符, 这也验证了所提出的基于系统性能衰退预测与黎曼度量的故障隔离方案的有效性.
值得注意的是, 当三容水箱处于不同工况时, 其性能预测矩阵也会不同. 因此, 可将不同的工况作为相应的性能变化或衰退模态考虑到分类中. 通过训练不同工况下可能的性能预测矩阵的均值和半径, 进而借助所提出的故障定位方法可实现对不同工况和不同类型的性能变化或衰退情况的区分. 因此, 本文所提出的方法不仅适用于系统老化、故障、控制器参数变化等因素引发的性能衰退, 也可用于系统不同工况的判定.
5. 总结与展望
针对一类带有反馈控制环节的动态系统, 本文提出了一种基于黎曼度量与控制性能衰退预测的性能监测与诊断方法. 首先, 提出了一个性能衰退的预测指标, 并给出该性能指标的离线与在线计算方法; 其次, 基于黎曼度量提出了系统性能衰退程度的监测方法, 并基于随机算法给出了对应的黎曼均值与阈值的设定方法, 进而实现性能监测; 最后, 通过分析各类故障的数据构建各类故障模态性能库并设计对应的阈值, 进而实现故障的实时定位. 通过三容水箱系统仿真验证了所提出故障诊断方法的有效性. 所提出的基于黎曼度量的性能监测与诊断方法既可以实现系统性能衰退监测, 又可识别故障位置, 在处理反馈控制系统的性能监测中显示出优秀的监测性能.
本文提出的性能监测与诊断方法面向线性系统. 但在该设计框架内, 可借助机器学习技术实现复杂工业非线性系统的性能衰退预测与故障诊断. 同时, 近年来基于性能的容错控制方法也受到越来越多的关注, 如何基于性能衰退预测的指标实现系统性能修复值得进一步研究. 另外, 目前工业系统结构复杂, 多呈现多子系统互联耦合的形式[44], 如何实现分布式性能衰退监测与容错控制, 也是值得研究的课题.
-
图 2 基于黎曼度量的控制性能监测流程图
Fig. 2 Flow chart of Riemannian metric based control performance monitoring
表 1 水箱DTS200的参数
Table 1 Parameters of tank DTS200
参数 符号 值 单位 水箱面积 $ {\cal{A}} $ 154 $ \mathrm{cm}^{2} $ 水箱间管道面积 $s_{{n} }$ 0.5 $ \mathrm{cm}^{2} $ 水箱最高水位 $ H_{\mathrm{max}} $ 62 $ \mathrm{cm} $ 泵 1 的最大进水量 $ Q_{1_{\mathrm{max}}} $ 120 $ \mathrm{cm}^{3}/\mathrm{s} $ 泵 2 的最大进水量 $ Q_{2_{\mathrm{max}}} $ 120 $ \mathrm{cm}^{3}/\mathrm{s} $ 管道 1 水流系数 $ a_{1} $ 0.45 管道 2 水流系数 $ a_{2} $ 0.60 管道 3 水流系数 $ a_{3} $ 0.45 
下载: 导出CSV
表 2 水箱堵塞故障隔离
Table 2 Isolation of pipe plugging
故障 $ d_R^2({\boldsymbol P},{\boldsymbol P}_{z,1}) $ $ d_R^2({\boldsymbol P},{\boldsymbol P}_{z,2}) $ $ d_R^2({\boldsymbol P},{\boldsymbol P}_{z,3}) $ 故障隔离 $a_1 = 0.30$ $ 0.0388 $ $ 0.1088 $ $ 0.1193 $ 故障簇 1 $ a_2 = 0.35 $ $ 0.0894 $ $ 0.0270 $ $ 0.0810 $ 故障簇 2 $ a_3 = 0.28 $ $ 0.1258 $ $ 0.1099 $ $ 0.0478 $ 故障簇 3 $ a_1 = 0.27 $ $ 0.0636 $ $ 0.1325 $ $ 0.1409 $ 故障簇 1 $a_2 = 0.40$ $ 0.0809 $ $ 0.0139 $ $ 0.0732 $ 故障簇 2
下载: 导出CSV
-
[1] 何潇, 郭亚琦, 张召, 贾繁林, 周东华. 动态系统的主动故障诊断技术. 自动化学报, 2020, 46(8): 1557-1570He Xiao, Guo Ya-Qi, Zhang Zhao, Jia Fan-Lin, Zhou Dong-Hua. Active fault diagnosis for dynamic systems. Acta Automatica Sinica, 2020, 46(8): 1557-1570 [2] Ding S X. Data-driven Design of Fault Diagnosis and Fault-tolerant Control Systems. Berlin: Springer, 2014. [3] 彭开香, 马亮, 张凯. 复杂工业过程质量相关的故障检测与诊断技术综述. 自动化学报, 2017, 43(3): 349-365Peng Kai-Xiang, Ma Liang, Zhang Kai. Review of quality-related fault detection and diagnosis techniques for complex industrial processes. Acta Automatica Sinica, 2017, 43(3): 349-365 [4] 钱峰, 杜文莉, 钟伟民, 唐漾. 石油和化工行业智能优化制造若干问题及挑战. 自动化学报, 2017, 43(6): 893-901Qian Feng, Du Wen-Li, Zhong Wei-Min, Tang Yang. Problems and challenges of smart optimization manufacturing in petrochemical industries. Acta Automatica Sinica, 2017, 43(6): 893-901 [5] Li L L, Ding S X. Gap metric techniques and their application to fault detection performance analysis and fault isolation schemes. Automatica, 2020, 118: Article No. 109029 [6] 陈晓露, 王瑞璇, 王晶, 周靖林. 基于混合型判别分析的工业过程监控及故障诊断. 自动化学报, 2020, 46(8): 1600-1614Chen Xiao-Lu, Wang Rui-Xuan, Wang Jing, Zhou Jing-Lin. Industrial process monitoring and fault diagnosis based on hybrid discriminant analysis. Acta Automatica Sinica, 2020, 46(8): 1600-1614 [7] Huang B, Kadali R. Dynamic modeling, predictive control and performance monitoring: A data-driven subspace approach. Lecture Notes in Control and Information Sciences. London: Springer, 2008. [8] Jelali M. An overview of control perormance assessment technology and industrial applications. Control Engineering Practice, 2006, 14: 441-466 doi: 10.1016/j.conengprac.2005.11.005 [9] Schafer J, Cinar A. Multivariable MPC system performance assessment, monitoring, and diagnosis. Journal of Process Control, 2004, 14(2): 113-129 doi: 10.1016/j.jprocont.2003.07.003 [10] Verenich I, Dumas M, Rosa M L, Nguyen H. Predicting process performance: A white-box approach based on process models. Journal of Software: Evolution and Process, 2019, 31(6): Article No. e2170 [11] Ding S X, Yin S, Peng K, Hao H. A novel scheme for key performance indicator prediction and diagnosis with application to an industrial hot strip mill. IEEE Transactions on Industrial Informatics, 2013, 9(4): 2239-2247 doi: 10.1109/TII.2012.2214394 [12] Xie X, Sun W, Cheung K. An advanced PLS approach for key performance indicator related prediction and diagnosis in case of outliers. IEEE Transactions on Industrial Electronics, 2016, 63(4): 2587-2594 [13] Duan Y, Liu M, Dong M. A metric-learning-based nonlinear modeling algorithm and its application in key-performance-indicator prediction. IEEE Transactions on Industrial Electronics, 2020, 67(8): 7073-7082 doi: 10.1109/TIE.2019.2935979 [14] Li L, Luo H, Ding S X, Yang Y, Peng K. Performance-based fault detection and fault-tolerant control for automatic control systems. Automatica, 2019, 99: 309-316 [15] Tao X, Lu C, Lu C, Wang Z. An approach to performance assessment and fault diagnosis for rotating machinery equipment. Eurasip Journal on Advances in Signal Processing, 2013, 2013(1): 1-16 doi: 10.1186/1687-6180-2013-1 [16] Li L, Ding S X. Performance supervised fault detection schemes for industrial feedback control systems and their data-driven implementation. IEEE Transactions on Industrial Informatics, 2020, 16(4): 2849-2858 doi: 10.1109/TII.2019.2940099 [17] Seem J E. Method and System for Assessing Performance of Control Systems, U.S. US7729882 B2, 2009. [18] Yin S, Zhu X, Kaynak O. Improved PLS focused on key-performance-indicator-related fault diagnosis. IEEE Transactions on Industrial Electronics, 2015, 62(3): 1651-1658 doi: 10.1109/TIE.2014.2345331 [19] He S, Wang Y, Liu C. Modified partial least square for diagnosing key-performance-indicator-related faults. The Canadian Journal of Chemical Engineering, 2018, 96(2): 444-454 doi: 10.1002/cjce.23002 [20] Song B, Zhou X, Shi H, Tao Y. Performance-indicator-oriented concurrent subspace process monitoring method. IEEE Transactions on Industrial Electronics, 2019, 66(7): 5535-5545 doi: 10.1109/TIE.2018.2868316 [21] Li H, Zhao J, Zhang X, Teng H. Gear fault diagnosis and damage level identification based on Hilbert transform and Euclidean distance technique. Journal of Vibroengineering, 2014, 16(8): 4137-4151 [22] Tian Y, Lu C, Wang Z. Self-adaptive bearing fault diagnosis based on permutation entropy and manifold-based dynamic time warping. Mechanical Systems & Signal Processing, 2019, 114: 658-673 [23] Patel S, Upadhyay S H. Euclidean distance based feature ranking and subset selection for bearing fault diagnosis. Expert Systems with Applications, 2020, DOI: 10.1016/j.eswa.2020.113400 [24] Liang H. Fault analysis of hierarchical cluster and fault diagnosis of Mahalanobis distance in analog circuit. Journal of Electronic Measurement and Instrument, 2010, 24(7): 610-615 doi: 10.3724/SP.J.1187.2010.00610 [25] 艾延廷, 冯研研, 周海仑.小波变换和EEMD-马氏距离的轴承故障诊断.噪声与振动控制, 2015, 1: 235-239 doi: 10.3969/j.issn.1006-1355.2015.01.050Ai Yan-Ting, Feng Yan-Yan, Zhou Hai-Lun. Fault diagnosis of roller bearings using wavelet transform and EEMD-Mahalanobis distance. Noise and Vibration Control, 2015, 1: 235-239\\ doi: 10.3969/j.issn.1006-1355.2015.01.050 [26] Ji H. Statistics Mahalanobis distance for incipient sensor fault detection and diagnosis. Chemical Engineering Sience, 2021, DOI: 10.1016/j.ces.2020.116233 [27] 吕鹏飞, 闫云聚, 荔越. 基于马氏距离的改进核Fisher化工故障诊断研究. 自动化学报, 2020, 46(11): 2379-2392Lv Peng-Fei, Yan Yun-Ju, Li Yue. Research on fault diagnosis of improved kernel Fisher based on Mahalanobis distance in the field of chemical industry. Acta Automatica Sinica, 2020, 46(11): 2379-2391 [28] Amari S. Information Geometry and Its Applications. Berlin: Springer, 2016. [29] Boothby W M. An Introduction to Differentiable Manifolds and Riemannian Geometry. Pittsburgh: Academic Press, 1975. [30] An J, Ai P. Deep domain adaptation model for bearing fault diagnosis with Riemann metric correlation alignment. Mathematical Problems in Engineering, 2020, 1: 1-12 [31] 周美含, 姜宏, 孙帅. 基于黎曼流形的MIMO雷达目标检测方法. 吉林大学学报: 信息科学版, 2020, 3: 237-242Zhou Mei-Han, Jiang Shuai, Suan Shuai. Target detection method for MIMO radar based on Riemannian manifold. Journal of Jilin University(Information Science Edition), 2020, 3: 237-242 [32] Wang S, Sun X, Li C. Wind turbine gearbox fault diagnosis method based on Riemannian manifold. Mathematical Problems in Engineering: Theory, Methods and Applications, 2014, 8: 1-10 [33] Wang Z, Jia L, Qin Y. Adaptive diagnosis for rotating machineries using information geometrical Kernel-ELM based on VMD-SVD. Entropy, 2018, 20(1): 73-91 doi: 10.3390/e20010073 [34] 孙小婷. 主测地线分析技术在汽轮机系统中的应用. 自动化应用, 2020, 1: 22-25Sun Xiao-Ting. Application of main geodesic analysis technology in steam turbine system. Automation Application, 2020, 1: 22-25 [35] Hiriart-Urruty J B, Malick J. A fresh variational-analysis look at the positive semidefinite matrices world. Journal of Optimization Theory and Applications, 2012, 153(3): 551-577 doi: 10.1007/s10957-011-9980-6 [36] Moakher M. A differential geometric approach to the geometric mean of symmetric positive-definite matrices. Siam Journal on Matrix Analysis & Applications, 2005, 26: 735-747 [37] Moakher M, Bathelor P G. Symmetric Positive-definite Matrices: From Geometry to Applications and Visualization. Berlin: Springer, 2006. [38] Ding S X. Advanced Methods for Fault Diagnosis and Fault-tolerant Control. Berlin: Springer, 2021. [39] Magnus J R. Linear Structures. Oxford: Oxford University Press, 1998. [40] Tempo R, Calcfiro G, Dabbene F. Randomized Algorithms for Analysis and Control of Uncertain Systems. Berlin: Springer, 2005. [41] Cryan M. Probability and computing: Randomized algorithms and probabilistic analysis. Bulletin of Symbolic Logic, 2006, 12(2): 304-308 doi: 10.1017/S107989860000278X [42] Ding S X, Li L, Krüger M. Application of randomized algorithms to assessment and design of observer-based fault detection systems. Automatica, 2019, 107: 175-182 doi: 10.1016/j.automatica.2019.05.037 [43] Lewis F L, Vrabie D, Vamvoudakis K G. Reinforcement learning and feedback control using natural decision methods to design optimal adaptive controllers. IEEE Control Systems Magazine, 2012, 32(6): 76-105 doi: 10.1109/MCS.2012.2214134 [44] 杨浩, 姜斌, 周东华. 互联系统容错控制的研究回顾与展望. 自动化学报, 2017, 43(1): 9-19Yang Hao, Jiang Bin, Zhou Dong-Hua. Review and perspectives on fault tolerant control for interconnected systems. Acta Automatica Sinica, 2017, 43(1): 9-19 期刊类型引用(0)
其他类型引用(2)
-
下载:
计量
- 文章访问数: 2555
- HTML全文浏览量: 499
- PDF下载量: 197
- 被引次数: 2
