基于扩展卡尔曼滤波的故障检测

打开实时脚本

这个例子展示了如何使用扩展卡尔曼滤波器进行故障检测。该实例使用扩展卡尔曼滤波器在线估计简单直流电机的摩擦。在估计的摩擦中检测到显著的变化，并表明一个故障。本示例使用系统识别工具箱™中的功能，不需要预测性维护工具箱™。

汽车模型

电机被建模为一个带有阻尼系数c的惯性J，由扭矩u驱动。电机的角速度w和加速度 $w_{}^{˙}$ ，是测量的输出。

$w_{}^{˙} ＝（ u - c w ） / J$

要使用扩展卡尔曼滤波器估计阻尼系数c，请为阻尼系数引入辅助状态，并将其导数设置为零。

$c_{}^{˙} ＝ 0$

因此，模型状态x = [w;c]和测量y的方程为:

$［ \begin{array}{c} w_{}^{˙} \\ c_{}^{˙} \end{array} ］＝［ \begin{array}{c} （ u - c w ） / J \\ 0 \end{array} ］$

$y ＝［ \begin{array}{c} w \\ （ u - c w ） / J \end{array} ］$

利用近似法将连续时间方程转化为离散时间方程 $x_{}^{˙} ＝ \frac{x_{n + 1} - x_{n}}{T_{年代}}$ ，其中Ts是离散采样周期。这给出了在pdmMotorModelStateFcn.m和pdmmotormodemeasurementfcn.m功能。

$［ \begin{array}{c} w_{n + 1} \\ c_{n + 1} \end{array} ］＝［ \begin{array}{c} w_{n} + （ u_{n} - c_{n} w_{n} ） T_{年代} / J \\ c_{n} \end{array} ］$

$y_{n} ＝［ \begin{array}{c} w_{n} \\ （ u_{n} - c_{n} w_{n} ） / J \end{array} ］$

指定电机参数。

J = 10;%的惯性Ts=0.01；%样品时间

指定初始状态。

x0 = [．．．0;．．．%角速度1);%的摩擦类型pdmMotorModelStateFcn

function x1 = pdmMotorModelStateFcn(x0,varargin) % pdmMotorModelStateFcn % %状态更新方程，以摩擦作为状态% % x1 = pdmMotorModelStateFcn(x0,u,J,Ts) % %% u -电机扭矩输入% J -电机惯性% Ts -采样时间% % output: % x1 - updated states % %% Input J = varargin{2};% System innertia Ts = varargin{3};% State update equation x1 =[…]x0 (1) + Ts / J * (u-x0 (1) * x0 (2));．．．x0 (2)];结束

类型pdmMotorModelMeasurementFcn

函数y=Pdmotromodelmeasurementfcn（x，varargin）%Pdmotromodelmeasurementfcn%%以摩擦为状态的电机的测量方程%%y=Pdmotromodelmeasurementfcn（x0，u，J，Ts）%%输入：%x-带元件的电机状态[角速度；摩擦]%u-电机扭矩输入%J-电机惯性%Ts-采样时间%%输出：%y-电机测量元件[角速度；角加速度]%2016-2017版权所有MathWorks，Inc.%从输入中提取数据u=varargin{1}；%input J=varargin{2}；%System innertia%输出方程y=[…x（1）；…（u-x（1）*x（2））/J]；结束

电机经历状态（过程）噪声干扰q和测量噪声干扰r。噪声项是相加的。

$［ \begin{array}{c} w_{n + 1} \\ c_{n + 1} \end{array} ］＝［ \begin{array}{c} w_{n} + （ u_{n} - c_{n} w_{n} ） T_{年代} / J \\ c_{n} \end{array} ］ + 问$

$y_{n} ＝［ \begin{array}{c} w_{n} \\ （ u_{n} - c_{n} w_{n} ） / J \end{array} ］ + r$

过程噪声和测量噪声均为零，E[问]= E [r] = 0，以及协方差Q = E (qq)和R=E[rr']．摩擦状态具有较高的过程噪声干扰。这反映了我们期望摩擦在电机正常运行时变化，并希望过滤器跟踪这种变化的事实。加速度和速度状态噪声很低，但速度和加速度测量相对有噪声。

指定过程噪声协方差。

Q = [．．．1 e-6 0;．．．%角速度0 1依照];%的摩擦

指定测量噪声协方差。

R = (．．．1 0的军医;．．．%速度测量0 1e-4]；%加速度测量

创建扩展卡尔曼滤波器

建立一个扩展的卡尔曼滤波来估计模型的状态。我们对阻尼状态特别感兴趣，因为这个状态值的剧烈变化表明故障事件。

创建一个扩展卡尔曼滤波器对象，并指定状态转移和测量函数的雅可比矩阵。

卡尔曼滤波器= extendedKalmanFilter (．．．@pdmMotorModelStateFcn，．．．@pdmMotorModelMeasurementFcn,．．．x0,．．．“状态协方差”, [1 0; 0 1000],．．．[1 0 0;0 1 0;00 100]，…“ProcessNoise”问,．．．“测量噪音”R．．．“StateTransitionJacobianFcn”@pdmMotorModelStateJacobianFcn,．．．“MeasurementJacobianFcn”, @pdmMotorModelMeasJacobianFcn);

扩展卡尔曼滤波器的输入参数是前面定义的状态转移和测量函数。初始状态值x0，初始状态协方差以及过程和测量噪声协方差也是扩展卡尔曼滤波器的输入。在本例中，可以从状态转移函数导出精确的雅可比函数f，测量功能h：

$\frac{\partial}{\partial x} f ＝［ \begin{array}{cc} 1 - T_{年代} c_{n} / J & - T_{年代} w_{n} / J \\ 0 & 1 \end{array} ］$

$\frac{\partial}{\partial x} h ＝［ \begin{array}{cc} 1 & 0 \\ - c_{n} / J & - w_{n} / J \end{array} ］$

状态雅可比矩阵定义在pdmMotorModelStateJacobianFcn.m函数和测量雅可比矩阵定义在pdmMotorModelMeasJacobianFcn.m函数。

类型pdmMotorModelStateJacobianFcn

function Jac = pdmMotorModelStateJacobianFcn(x,varargin) % pdmMotorModelStateJacobianFcn % %电机模型状态方程的雅可比矩阵。请参阅pdmMotorModelStateFcn获取%模型方程。% % Jac = pdmMotorModelJacobianFcn(x,u,J,Ts) % % Inputs: % x - state with elements[角速度;% u -电机扭矩输入% J -电机惯性% Ts -采样时间% % output: % jacc - state Jacobian computed at x % %Ts =变长度输入宗量{3};% Jacobian Jac =[…]1 - t / J * (2) ts / J * x (1);．．．0 1];结束

类型pdmMotorModelMeasJacobianFcn

函数J = pdmMotorModelMeasJacobianFcn(x,varargin) % pdmMotorModelMeasJacobianFcn % %电机模型测量方程的雅可比矩阵。参见% pdmMotorModelMeasurementFcn获取模型方程。% % Jac = pdmMotorModelMeasJacobianFcn(x,u,J,Ts) % % Inputs: % x - state with elements[角速度;% u -电机扭矩输入% J -电机惯性% Ts -采样时间% % output: % jacc -测量雅可比矩阵计算在x % %% System innertia % Jacobian J =[…]1 0;J - x - x (2) / (1) / J];结束

模拟

为了模拟工厂，创建一个回路，并在电机中引入一个故障(在电机小说中是一个戏剧性的变化)。在仿真回路中，使用扩展卡尔曼滤波器估计电机状态，并具体跟踪摩擦状态，检测何时有统计上显著的摩擦变化。

该电机是模拟脉冲序列，重复加速和减速的电机。这种类型的电机操作是生产线上拣选机器人的典型操作。

t = 0: Ts: 20;%时间，20秒，带Ts采样周期u =双(mod (t, 1) < 0.5) -0.5;%脉冲序列，周期1,50%占空比元=元素个数(t);%时间点个数nx =大小(x0, 1);%状态数ySig = 0 ([2, nt]);被测电机输出百分比xSigTrue = 0 ([nx, nt]);%未测量的运动状态xSigEst = 0 ([nx, nt]);运动状态估计XSTD = 0 ([nx nx nt]);%估计状态的标准差ySigEst = zeros([2, nt]);模型估计输出fMean = 0(1元);%平均估计摩擦力fSTD = 0(1元);%估计摩擦的标准偏差fKur = 0(2元);%估计摩擦峰度fChanged=false（1，nt）；%指示摩擦变化检测的标志

在模拟电机时，添加与构建扩展卡尔曼滤波器时使用的Q和R噪声协方差值类似的过程和测量噪声。对于摩擦，使用更小的噪声值，因为除了故障发生时，摩擦基本上是恒定的。在模拟过程中人为诱发故障。

rng (“默认”)；Qv=chol（Q）；%过程噪声的标准偏差Qv(结束)= 1飞行;%摩擦噪声小Rv=chol（R）；%测量噪音的标准偏差

使用状态更新方程模拟模型，并在模型状态中添加过程噪声。进入模拟十秒后，迫使马达摩擦力发生变化。使用模型测量功能模拟电机传感器，并在模型输出中添加测量噪声。

为ct=1:numel（t）型号输出更新y=PdM模型测量Cn（x0，u（ct），J，Ts）；y=y+Rv*randn（2,1）；添加测量噪声ySig（：，ct）=y；模型状态更新xSigTrue (:, ct) = x0;x1 = pdmMotorModelStateFcn (x0, u (ct), J, Ts);引起摩擦变化如果T (ct) == 10 x1(2) = 10;%阶跃变化结束x1n=x1+Qv*randn（nx，1）；%添加过程噪声x1n（2）=最大值（x1n（2），0.1）；摩擦下限x0 = x1n;%存储下一次模拟迭代的状态

从运动测量来估计运动状态，使用预测和正确的扩展卡尔曼滤波器的命令。

使用扩展卡尔曼滤波器的状态估计x_corr =正确(卡尔曼滤波器,y, u (ct), J, Ts);修正基于当前测量的状态估计。xSigEst (:, ct) = x_corr;xstd (:,:, ct) =胆固醇(ekf.StateCovariance);预测(ekf, u (ct), J, Ts);根据当前状态和输入预测下一个状态。

为了检测摩擦的变化，使用4秒移动窗口计算估计摩擦平均值和标准偏差。在最初的7秒周期后，锁定计算的平均值和标准偏差。这个初始计算的平均值是期望的无故障摩擦平均值。7秒后，如果估计的摩擦力距离预期的无故障平均值大于3个标准差，则表示摩擦力发生了显著变化。为了减少噪声和估计摩擦力变异性的影响，在与3个标准偏差界限比较时使用估计摩擦力的平均值。

如果t (ct) < 7%计算估计值的平均值和标准偏差。idx = max (ct - 400):马克斯(1,ct-1);% Ts = 0.01秒fMean（ct）=平均值（xSigEst（2，idx））；fSTD（ct）=标准值（xSigEst（2，idx））；其他的%存储计算的平均值和标准偏差，无需%验算。fMean (ct) = fMean (ct-1);fSTD (ct) = fSTD (ct-1);%用期望的摩擦平均值和标准偏差来检测%摩擦的变化。estFriction =意味着(xSigEst(2,马克斯(1,ct-10): ct));fchange (ct) = (estFriction > fMean(ct)+3*fSTD(ct)) || (estFriction < fMean(ct)-3*fSTD(ct));结束如果fChanged（ct）和&~fChanged（ct-1）%在摩擦变化信号|fChanged|中检测到上升边。流('在%f\n处有显著的摩擦变化'，t（ct））；结束

在10.450000处有显著的摩擦变化

使用估计状态来计算估计输出。计算测量输出和估计输出之间的误差，并计算误差统计量。误差统计量可用于检测摩擦变化。稍后将对此进行更详细的讨论。

ySigEst (:, ct) = pdmMotorModelMeasurementFcn (x_corr u (ct), J, Ts);idx = max (ct - 400): ct;fKur (:, ct) = (．．．idx峰度(ySigEst (1) -ySig (idx));．．．峰度（ySigEst（2，idx）-ySig（2，idx））；结束

扩展卡尔曼滤波性能

请注意，在10.45秒时检测到摩擦变化。我们现在描述如何推导此故障检测规则。首先检查模拟结果和过滤器性能。

图，子地块（211），地块（t，ySig（1，：），t，ySig（2，：）；头衔(“电机输出”)传说(测量角速度的，测量角加速度的，“位置”，“西南”)次要情节(212)、情节(t, u);标题(“电机输入-扭矩”）

图中包含2个轴对象。具有标题电机输出的轴对象1包含2个line类型的对象。这些对象表示测量的角速度、测量的角加速度。具有标题电机输入-扭矩的轴对象2包含line类型的对象。

模型的输入输出响应表明，直接从测量信号中检测摩擦变化是困难的。扩展卡尔曼滤波使我们能够估计状态，特别是摩擦状态。比较真实模型状态和估计状态。估计状态用3个标准差对应的置信区间表示。

图中,次要情节(211),情节(t, xSigTrue (1:), t, xSigEst (1:)．．．[t南t]、[xSigEst(1:) + 3 *挤压(xstd(1 1:))”,南xSigEst(: 1) 3 *挤压(xstd(1 1:))的])轴(20 -0.06 - 0.06[0]),传说(“真实价值”，“估计价值”，置信区间的)标题(“电机状态-速度”)子地块（212），地块（t，xSigTrue（2，：），t，xSigEst（2，：），．．．[t南t]、[xSigEst(2:) + 3 *挤压(xstd(2, 2,:))的南xSigEst(2:) 3 *挤压(xstd(2, 2,:))的])轴([0 -10 15])标题(“运动状态-摩擦”)；

图中包含2个轴对象。标题为“电机状态-速度”的轴对象1包含3个类型为line的对象。这些对象分别表示真值、估计值、置信区间。轴对象2，标题为电机状态-摩擦包含3个类型为线的对象。

注意，过滤器估计跟踪真值，并且置信区间保持有界。检查估计错误可以更深入地了解过滤器的行为。

图中,次要情节(211),情节(t, xSigTrue (1:) -xSigEst(1:))标题(“速度状态错误”)子地块（212），地块（t，xSigTrue（2，：）-xSigEst（2，：））标题(“摩擦状态错误”）

图中包含2个轴对象。标题为“速度状态错误”的轴对象1包含一个类型为line的对象。标题为“摩擦状态错误”的轴对象2包含一个类型为line的对象。

误差图表明，该滤波器在摩擦变化10秒后自适应，并将估计误差降至零。然而，错误图不能用于故障检测，因为它们依赖于知道真实状态。将测量到的状态值与加速度和速度的估计状态值进行比较可以提供一种检测机制。

图子地块（211），地块（t，ySig（1，：）-ySigEst（1，：）标题(“速度测量误差”)次要情节(212)、图(t, ySig (2:) -ySigEst(2:))标题(加速度测量误差的）

图中包含2个轴对象。轴对象1标题为速度测量误差包含一个类型线对象。标题为“加速度测量错误”的轴对象2包含一个类型为line的对象。

加速度误差图显示，当故障引入时，平均误差在10秒左右有微小差异。查看错误统计信息，看看是否可以从计算出的错误中检测到故障。加速度和速度误差被期望是正态分布(噪声模型都是高斯)。因此，加速度误差的峰度可以帮助识别摩擦力变化引起的误差分布由对称向非对称变化，从而引起误差分布的变化。

图中,次要情节(211),情节(t, fKur(1:))标题(“速度误差峰度”)次要情节(212)、图(t, fKur(2:))标题(“加速度误差峰度”）

图中包含两个轴对象。具有标题速度误差峰度的轴对象1包含类型为line的对象。具有标题加速度错误峰度的轴对象2包含类型为line的对象。

忽略前4秒，当估计器仍在收敛和数据收集，误差的峰度是相对恒定的，微小的变化在3(高斯分布的期望峰度值)。因此，在这个应用程序中，错误统计不能用于自动检测摩擦变化。在这个应用程序中，使用误差的峰度也是困难的，因为滤波器是自适应的，并不断地将误差驱动到零，只给出一个误差分布不为零的短时间窗口。

因此，在这种应用中，使用估计摩擦的变化提供了自动检测电机故障的最佳方法。来自已知无故障数据的摩擦估计(平均值和标准偏差)提供了摩擦的预期边界，当这些边界被违反时，很容易检测到。下图强调了这种故障检测方法。

图(t,xSigEst(2，:)，[t nan t]，[fMean+3*fSTD,nan,fMean-3*fSTD])摩擦变化检测的)传说(“估计摩擦”，“无过错摩擦边界”)轴([0 20 -10 20])网格在…上

图中包含一个轴对象。标题为“摩擦变化检测”的轴对象包含2个line类型的对象。这些对象表示估计的摩擦，无故障摩擦边界。

总结

这个例子展示了如何使用扩展卡尔曼滤波器估计摩擦在一个简单的直流电动机和使用摩擦估计故障检测。

基于扩展卡尔曼滤波的故障检测

汽车模型

创建扩展卡尔曼滤波器

模拟

扩展卡尔曼滤波性能

总结

相关话题

预测性维护工具箱文档

金宝app

基于MATLAB的预测维护