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

打开实时脚本

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

汽车模型

电机被建模为一个惯量J和阻尼系数c，由力矩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而且pdmMotorModelMeasurementFcn.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) % %输入:% x0 -初始状态与元素[角速度;摩擦]% u -电机转矩输入% J -电机惯量% Ts -采样时间% %输出:% x1 -更新状态% %版权所有2016-2017 The MathWorks, Inc. %从输入中提取数据u = varargin{1};%输入J = varargin{2};%系统内部Ts = varargin{3};状态更新方程x1 =[…]x0 (1) + Ts / J * (u-x0 (1) * x0 (2));．.．x0 (2)];结束

类型pdmMotorModelMeasurementFcn

function y = pdmMotorModelMeasurementFcn(x,varargin) % pdmMotorModelMeasurementFcn % %以摩擦为状态的电机测量方程% % y = pdmMotorModelMeasurementFcn(x0,u,J,Ts) % %输入:% x -电机状态与元素[角速度;% u -电机转矩输入% J -电机惯量% Ts -采样时间% %输出:% y -电机测量与元件[角速度;The MathWorks, Inc. %从输入中提取数据u = varargin{1};%输入J = varargin{2};输出方程y =[…]x (1);．.．(u-x (1) * (2) x) / 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 1的军医];加速度测量

创建一个扩展卡尔曼滤波器

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

创建一个extendedKalmanFilter对象，并指定状态转换函数和测量函数的雅可比矩阵。

ekf = extendedKalmanFilter(.．.@pdmMotorModelStateFcn,.．.@pdmMotorModelMeasurementFcn,.．.x0,.．.“StateCovariance”， [10 0;0 1000),.．.[1 0 0;0 10 0;00 100]…“ProcessNoise”问,.．.“MeasurementNoise”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) % %输入:% x -状态与元素[角速度;% u -电机转矩输入% J -电机惯量% Ts -采样时间% %输出:% Jac -状态雅可比矩阵计算在x % %版权所有2016-2017 The MathWorks, Inc. %模型属性J = varargin{2};Ts = varargin{3};雅可比矩阵=[…]1 - t / J * (2) ts / J * x (1);．.．0 1];结束

类型pdmMotorModelMeasJacobianFcn

function J = pdmMotorModelMeasJacobianFcn(x,varargin) % pdmMotorModelMeasJacobianFcn % %电机模型测量方程的雅可比矩阵。模型方程见% pdmMotorModelMeasurementFcn。% % Jac = pdmMotorModelMeasJacobianFcn(x,u,J,Ts) % %输入:% x -状态与元素[角速度;% u -电机转矩输入% J -电机惯量% Ts -采样时间% %输出:% Jac -测量雅可比矩阵计算在x % %版权所有2016-2017 The MathWorks, Inc. %系统参数J = varargin{2};%系统内部%雅可比矩阵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%Nt =数字(t);%时间点个数Nx = size(x0,1);%状态数ySig = 0 ([2, nt]);测量电机输出xSigTrue = 0 ([nx, nt]);%未测量的电机状态xSigEst = 0 ([nx, nt]);%估计运动状态XSTD = 0 ([nx nx nt]);%估计状态的标准差ySigEst = 0 ([2, nt]);%估计模型输出fMean = 0 (1,nt);平均估计摩擦fSTD = 0 (1,nt);估计摩擦的标准偏差fKur = 0 (2,nt);估计摩擦的峰度fChanged = false(1,nt);%表示摩擦变化检测的标志

在模拟电机时，添加与构建扩展卡尔曼滤波器时使用的Q和R噪声协方差值相似的过程噪声和测量噪声。对于摩擦，使用一个小得多的噪声值，因为摩擦在大多数情况下是恒定的，除非发生故障。模拟过程中人为诱发故障。

rng (“默认”）;Qv = chol(Q);过程噪声的标准偏差Qv(end) = 1e-2;%摩擦噪声更小Rv = chol(R);测量噪声的标准偏差

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

为Ct = 1:数字(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 = correct(ekf,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(1,ct-400):max(1,ct-1);% Ts = 0.01秒fMean(ct) = mean(xSigEst(2, idx));fSTD(ct) = std(xSigEst(2, idx));其他的存储计算的平均值和标准偏差%验算。fMean(ct) = fMean(ct-1);fSTD(ct) = fSTD(ct-1);使用期望摩擦平均值和标准偏差进行检测%摩擦变化。estFriction = mean(xSigEst(2,max(1,ct-10):ct));fChanged(ct) = (estFriction > fMean(ct)+3*fSTD(ct)) || (estFriction < fMean(ct)-3*fSTD(ct));结束如果fChanged(ct) && ~fChanged(ct-1)在摩擦变化信号|fChanged|中检测上升边。流('在%f\n处的显著摩擦变化't (ct));结束

在10.45万处有显著摩擦变化

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

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

扩展卡尔曼滤波性能

注意，在10.45秒时检测到摩擦变化。现在我们将描述这个故障检测规则是如何导出的。首先检验仿真结果和滤波器性能。

图中,次要情节(211),情节(t, ySig (1:), t, ySig (2:));标题(电机输出的)传说(“测量角速度”，“测量角加速度”，“位置”，“西南”) subplot(212)， plot(t,u);标题(“电机输入-转矩”）

图中包含2个轴对象。axis对象1的标题为Motor Outputs，包含2个类型为line的对象。这些物体代表测得的角速度，测得的角加速度。标题为Motor Input - Torque的Axes对象2包含一个类型为line的对象。

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

图，subplot(211)，plot(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个轴对象。标题为Motor State - Velocity的轴对象1包含3个类型为line的对象。这些对象表示真实值、估计值、置信区间。轴对象2，标题为Motor State -摩擦力，包含3个类型为line的对象。

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

图，subplot(211)，plot(t,xSigTrue(1，:)-xSigEst(1，:))“速度状态误差”)次要情节(212)、图(t, xSigTrue (2:) -xSigEst(2:))标题(“摩擦状态误差”）

图中包含2个轴对象。标题为Velocity State Error的Axes对象1包含一个类型为line的对象。标题为摩阻状态错误的Axes对象2包含一个类型为line的对象。

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

图subplot(211)， plot(t,ySig(1，:)-ySigEst(1，:))“速度测量误差”)次要情节(212)、图(t, ySig (2:) -ySigEst(2:))标题(“加速度测量误差”）

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

加速度误差图显示，当故障出现时，10秒左右的平均误差略有不同。查看错误统计信息，以查看是否可以从计算错误中检测到故障。加速度和速度误差预期为正态分布(噪声模型均为高斯分布)。因此，加速度误差的峰度可以帮助识别由于摩擦力的变化导致误差分布从对称到不对称的变化。

图，subplot(211)，plot(t,fKur(1，:))“速度误差峰度”) subplot(212)，plot(t,fKur(2，:))加速度误差峰度）

图中包含2个轴对象。标题为Velocity Error Kurtosis的Axes对象1包含一个类型为line的对象。标题为Acceleration Error Kurtosis的Axes对象2包含一个类型为line的对象。

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

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

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

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

总结

本例展示了如何使用扩展卡尔曼滤波器估计简单直流电动机的摩擦，并将摩擦估计用于故障检测。

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

汽车模型

创建一个扩展卡尔曼滤波器

模拟

扩展卡尔曼滤波性能

总结

相关的话题