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

打开直播脚本

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

电机模型

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

${W.}_{}^{˙} = （你 - C W. ） / j$

为了使用扩展的卡尔曼滤波器来估计阻尼系数C，引入阻尼系数的辅助状态，并将其导数设定为零。

$C_{}^{˙} = 0.$

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

$[\begin{array}{c} {W.}_{}^{˙} \\ C_{}^{˙} \end{array}] = [\begin{array}{c} （你 - C W. ） / j \\ 0. \end{array}]$

$y = [\begin{array}{c} W. \\ （你 - C W. ） / j \end{array}]$

使用近似将连续时间方程转换为离散时间 $X_{}^{˙} = \frac{X_{N + 1} - X_{N}}{{T.}_{S.}}$ 式中，Ts为离散采样周期。给出了离散时间模型方程，并在pdmmotormodelstatefcn.m.和pdmmotormodelmeasurementfcn.m.功能。

$[\begin{array}{c} {W.}_{N + 1} \\ C_{N + 1} \end{array}] = [\begin{array}{c} {W.}_{N} + （你_{N} - C_{N} {W.}_{N} ） {T.}_{S.} / j \\ C_{N} \end{array}]$

$y_{N} = [\begin{array}{c} {W.}_{N} \\ （你_{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 = pdmmotormodelmeasurementfcn（x，varargin）％pdmmotormodelmeasurementfcn％的电动机的测量方程，其具有摩擦的电动机作为状态％y = pdmmotormodelmeasurementfcn（x0，u，j，ts）％epce：％x  - 带元素的电机状态[角度速度;摩擦]％U  - 电动机扭矩输入％J  - 电动机惯量％TS  - 采样时间%%输出：％Y  - 电机测量与元素[角速度;角度加速度]％％2016-2017 MathWorks，Inc.％从输入中提取数据u = varargin {1};% Input J = varargin{2};％system interia％输出方程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} + （你_{N} - C_{N} {W.}_{N} ） {T.}_{S.} / j \\ C_{N} \end{array}] + 问：$

$y_{N} = [\begin{array}{c} {W.}_{N} \\ （你_{N} - C_{N} {W.}_{N} ） / j \end{array}] + R.$

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

指定过程噪声协方差。

q = [．..1E-6 0;．..％ 角速度0 1依照];%的摩擦

指定测量噪声协方差。

r = [．..1 0的军医;．..%速度测量0 1E-4];％加速度测量

创建扩展卡尔曼滤波器

创建一个扩展的卡尔曼筛选器以估计模型的状态。我们对阻尼状态特别感兴趣，因为此状态值的急剧变化表示故障事件。

创建一个ExtendedKalmanFilter.对象，并指定状态转换和测量功能的jacobians。

EKF = ExtendedKalmanFilter（．..@pdmmotormodelstatefcn，．..@pdmMotorModelMeasurementFcn,．..x0，．..'statecovariance'，[1 0;0 1000]，．..[1 0 0;0 1 0;00 100]，…'processnoise'问,．..'MeasurementNoise'R．..'stateTransitionJacobianfcn'@pdmMotorModelStateJacobianFcn,．..'measurementjacobianfcn', @pdmMotorModelMeasJacobianFcn);

扩展卡尔曼滤波器的输入参数是前面定义的状态转移和测量函数。初始状态值x0，初始状态协方差以及流程和测量噪声CovariRces也输入扩展的卡尔曼滤波器。在此示例中，可以从状态转换函数导出精确的Jacobian函数F，测量功能H：

$\frac{\partial}{\partial X} F = [\begin{array}{cc} 1 - {T.}_{S.} C_{N} / j & - {T.}_{S.} {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;％时间，20S有TS采样周期U = DOUBLE（MOD（T，1）<0.5）-0.5;%脉冲序列，周期1,50%占空比元=元素个数(t);%时间点个数nx = size（x0,1）;%状态数ySig = 0 ([2, nt]);被测电机输出百分比xSigTrue = 0 ([nx, nt]);％未测量的电机状态xSigEst = 0 ([nx, nt]);％估计的电机状态xstd = zeros（[nx nx nt]）;估计州的％标准偏差ysigest = zeros（[2，nt]）;模型估计输出fMean = 0(1元);％平均估计摩擦fstd = zeros（1，nt）;估计摩擦的％标准偏差FKUR =零（2，NT）;估计摩擦的％Kurtosis效率=假（1，NT）;％标志表示摩擦变化检测

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

rng (“默认”）;qv = chol（q）;工艺噪声的％标准偏差Qv(结束)= 1飞行;摩擦噪声较小rv = chol（r）;测量噪声的％标准偏差

使用状态更新方程式模拟模型，并将进程噪声添加到模型状态。在模拟中十秒钟，迫使电机摩擦的变化。使用模型测量功能模拟电机传感器，并向型号输出添加测量噪声。

为了ct = 1：numel（t）型号输出更新y = pdmmotormodelmeasurementfcn（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）= max（x1n（2），0.1）;摩擦下限x0 = x1n;%存储下一次模拟迭代的状态

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

使用扩展卡尔曼滤波器的状态估计x_corr =正确(卡尔曼滤波器,y, u (ct), J, Ts);％基于当前测量校正状态估计。xsigest（：，ct）= x_corr;xstd（：，：，ct）= chol（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）= STD（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));结尾如果效率（CT）&&〜效率（CT-1）％检测摩擦变化信号中的上升沿|效率。fprintf（'％f \ n'的显着摩擦变化，t（ct））;结尾

在10.450000处有显著的摩擦变化

使用估计状态计算估计的输出。计算测量和估计的输出之间的错误，并计算错误统计信息。错误统计数据可用于检测摩擦变化。稍后将更详细地讨论这一点。

ySigEst (:, ct) = pdmMotorModelMeasurementFcn (x_corr u (ct), J, Ts);idx = max (ct - 400): ct;fKur (:, ct) = (．..kurtosis（ysigest（1，Idx）-ysig（1，Idx））;．..kurtosis（ysigest（2，Idx） -  agsig（2，Idx））];结尾

扩展卡尔曼过滤性能

注意，在10.45秒内检测到摩擦变化。我们现在描述了如何派生此故障检测规则。首先检查仿真结果和过滤性能。

图，子图（211），plot（t，ysig（1，:)，t，ysig（2，:)）;标题（'电机输出')传说('测量角速度'那测量角加速度的那'地点'那“西南”）子图（212），绘图（t，u）;标题（'电机输入 - 扭矩'）

图中包含2个轴。带标题电机输出的轴1包含2个类型的线。这些对象表示测量的角速度，测得的角速度。具有标题电机输入的轴2 - 扭矩包含类型线的物体。

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

图中,次要情节(211),情节(t, xSigTrue (1:), t, xSigEst (1:)．..[t nan t]，[xsigest（1，:) + 3 *挤压（xstd（1,1，:)）'，nan，xsigest（1，:)  -  3 *挤压（xstd（1,1，:)）']）轴（[0 20 -0.06 0.06]），传奇（“真实价值”那“估计价值”那'置信区间')标题('电机状态 - 速度'）子平板（212），绘图（t，xsigtrue（2，:)，t，xsigest（2，:)，．..[t nan t]，[xsigest（2，:) + 3 *挤压（xstd（2,2，:)）'nan xsigest（2，:)  -  3 *挤压（xstd（2,2，:)'']）轴（[0 20 -10 15]）标题（'电机状态 - 摩擦'）;

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

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

图，子图（211），plot（t，xsigtrue（1，：） -  xsigest（1，:)）标题（“速度状态错误”）子图（212），plot（t，xsigtrue（2，：） -  xsigest（2，:)）标题（'摩擦状态错误'）

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

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

图形子图（211），plot（t，ysig（1，:)  -  ysigest（1，:)）标题（'速度测量错误')次要情节(212)、图(t, ySig (2:) -ySigEst(2:))标题(加速度测量误差的）

图中包含2个轴。具有标题速度测量误差的轴1包含类型线的对象。带标题加速度测量误差的轴2包含类型线的对象。

加速度误差绘图显示出故障时的2秒左右的误差差异。查看错误统计信息，以查看是否可以从计算的错误中检测到故障。预期加速度和速度误差通常是分布（噪声模型都是高斯）。因此，加速度误差的Kurtosis可以帮助识别误差分布由于摩擦变化而从对称变为不对称的变化，并且误差分布的变化。

图中,次要情节(211),情节(t, fKur(1:))标题(“速度误差峰度”）子图（212），绘图（t，fkur（2，:)）标题（'加速误差kurtosis'）

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

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

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

图(t,xSigEst(2，:)，[t nan t]，[fMean+3*fSTD,nan,fMean-3*fSTD])摩擦变化检测的)传说('估计摩擦'那'无故障摩擦界')轴([0 20 -10 20])网格在

图包含轴。具有标题摩擦变化检测的轴包含2个类型的类型。这些对象代表估计的摩擦，无故障摩擦界限。

总结

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

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

电机模型

创建扩展卡尔曼滤波器

模拟

扩展卡尔曼过滤性能

总结

相关话题

预测维护工具箱文档

万博1manbetx

基于MATLAB的预测维护