使用扩展卡尔曼滤波器进行故障检测

打开直播脚本

此示例显示如何使用扩展的卡尔曼滤波器进行故障检测。该示例使用扩展的卡尔曼滤波器进行在线估计简单直流电机的摩擦。检测估计摩擦的显着变化并指示故障。此示例使用系统识别工具箱™的功能，并且不需要预测性维护工具箱™。

汽车模型

电动机被建模为具有阻尼系数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.和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];％摩擦typepdmmotormodelstatefcn.

功能x1 = pdmmotormodelstatefcn（x0，varargin）％pdmmotormodelstatefcn％upsytey状态更新的电机与摩擦为一个状态%% x1 = pdmmotormodelstatefcn（x0，u，j，ts）％inclue：％x0  - 具有元素的初始状态[角速度;摩擦]％U  - 电动机扭矩输入％J  - 电动机惯量％TS  - 采样时间%%输出：％X1  - 更新状态％％2016-2017 MathWorks，Inc.％从输入中提取数据u = varargin {1};％输入j = varargin {2};％system innertia ts = varargin {3};％采样时间％状态更新等式x1 = [... x0（1）+ ts / j *（u-x0（1）* x0（2））;... x0（2）];结束

typepdmmotormodelmeasurementfcn.

功能y = pdmmotormodelmeasurementfcn（x，varargin）％pdmmotormodelmeasurementfcn％的电动机的测量方程，其具有摩擦的电动机作为状态％y = pdmmotormodelmeasurementfcn（x0，u，j，ts）％epce：％x  - 带元素的电机状态[角度速度;摩擦]％U  - 电动机扭矩输入％J  - 电动机惯量％TS  - 采样时间%%输出：％Y  - 电机测量与元素[角速度;角度加速度]％％2016-2017 MathWorks，Inc.％从输入中提取数据u = varargin {1};％输入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} + （ 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 [q] = e [r] = 0和coveramcesq = e [qq']和r = e [rr']．摩擦状态具有高过程噪声干扰。这反映了我们预期摩擦在电动机的正常操作期间变化的事实，并且希望过滤器跟踪这种变化。加速度和速度状态噪声低，但速度和加速度测量相对嘈杂。

指定过程噪声协方差。

Q = [......1 e-6 0;......%角速度0 1E-2];％摩擦

指定测量噪声协方差。

R = (......1E-4 0;......％速度测量0 1E-4];％加速度测量

创建扩展的卡尔曼过滤器

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

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

卡尔曼滤波器= extendedKalmanFilter (......@pdmmotormodelstatefcn，......@pdmmotormodelmeasurementfcn，......x0,......'statecovariance'，[1 0;0 1000]，......[1 0 0;0 1 0;0 0 100]，...“ProcessNoise”，q，......'MeasurementNoise'，r，......“StateTransitionJacobianFcn”，@pdmotormodelstatejacobianfcn，......“MeasurementJacobianFcn”，@pdmotormodelmeasjacobianfcn）;

扩展的卡尔曼滤波器具有输入参数，先前定义的状态转换和测量函数。初始状态值X0.，初始状态协方差以及流程和测量噪声CovariRces也输入扩展的卡尔曼滤波器。在此示例中，可以从状态转换函数导出精确的Jacobian函数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.函数。

typepdmmotormodelstatejacobianfcn.

功能Jac = PDMMotormodelStateJacobianFCN（x，varargin）％pdmmotormodelstatejacobianfcn％jacobian电机模型状态方程。请参阅pdmmotormodelstatefcn for％模型方程。%% jac = pdmmotormodeljacobianfcn（x，u，j，ts）%%输入：％x  - 具有元素[角速度;摩擦]％U  - 电动机扭矩输入％J  - 电动机惯量％TS  - 采样时间%%产出：％Jac  - 状态雅比尼亚在X％％的x％％2016-2017计算MathWorks，Inc.％模型属性J = varargin {2};ts = varargin {3};％jacobian jac = [... 1-ts / j * x（2）-ts / j * x（1）;... 0 1];结束

typepdmmotormodelmeasjacobianfcn.

功能J = PDMMOTORMODELMEAJACOBIANFCN（X，VARARGI）％PDMMOTORMODELMEASJACOBIANFCN％jacobian电机模型测量方程。有关模型方程，请参阅％pdmmotormodelmeasurementfcn。％jac = pdmmotormodelmeasjacobianfcn（x，u，j，ts）%%输入：％x  - 具有元素[角速度;摩擦]％U  - 电动机扭矩输入％J  - 电动机惯量％TS  - 采样时间％%%产出：％Jac  - 测量雅各比斯在X％％版权所有2016-2017 MathWorks，Inc.％系统参数J = varargin {2}。;％system interia％jacobian j = [... 1 0;-x（2）/ j -x（1）/ j];结束

模拟

为了模拟工厂，在电机（电动机小说中的剧烈变化）中产生循环并引入故障。在仿真循环中，使用扩展的卡尔曼滤波器来估计电机状态，并具体跟踪摩擦状态以检测摩擦统计上有显着变化。

用脉冲序列模拟电动机，该脉冲序列重复加速和减速电动机。这种类型的电机操作对于生产线中的选择器机器人来说是典型的。

T = 0：TS：20;％时间，20S有TS采样周期u =双(mod (t, 1) < 0.5) -0.5;％脉冲火车，第1期，50％占空比nt = numel（t）;％时间点数nx =大小(x0, 1);％的国家数量ysig = zeros（[2，nt]）;％测量电机输出xsigtrue = zeros（[nx，nt]）;%未测量的运动状态xsigest =零（[nx，nt]）;运动状态估计XSTD = 0 ([nx nx nt]);%估计状态的标准差ySigEst = zeros([2, nt]);％估计模型输出fmean =零（1，nt）;％平均估计摩擦fSTD = 0(1元);%估计摩擦的标准偏差fKur = 0(2元);估计摩擦的％Kurtosis效率=假（1，NT）;％标志表示摩擦变化检测

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

RNG（'默认'）;qv = chol（q）;工艺噪声的％标准偏差QV（END）= 1E-2;摩擦噪声较小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 =正确（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）=均值（xsigest（2，Idx））;FSTD（CT）= STD（Xsigest（2，IDX））;其他的％存储计算的均值和标准偏差％重新计算。FMEAN（CT）= FMEAN（CT-1）;FSTD（CT）= FSTD（CT-1）;％使用预期的摩擦均值和标准偏差来检测％摩擦变化。estfriction =平均值（xsigest（2，max（1，ct-10）：ct））;效率（CT）=（estfriction> fmean（ct）+ 3 * fstd（ct））||（estfriction 结束如果效率（CT）&&〜效率（CT-1）%在摩擦变化信号|fChanged|中检测到上升边。流('在%f\n处有显著的摩擦变化'，t（ct））;结束

10.450000的显着摩擦变化

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

ysigest（：，ct）= pdmmotormodelmeasurementfcn（x_corr，u（ct），j，ts）;IDX = MAX（1，CT-400）：CT;FKUR（：，CT）= [......idx峰度(ySigEst (1) -ySig (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），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个轴。具有标题电机状态 - 速度的轴1包含3个类型的线。这些对象表示真值，估计值，置信区间。具有标题电机状态 - 摩擦的轴2包含3个类型的线。

请注意，过滤器估计追踪真值，并且置信区间保持界限。检查估计错误提供更详细的过滤器行为。

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

图包含2个轴。具有标题速度状态误差的轴1包含类型线的对象。标题摩擦状态误差的轴2包含类型线的对象。

错误图显示过滤器在摩擦变化时适应10秒后，并将估计误差减少为零。但是，错误绘图不能用于故障检测，因为它们依赖于了解真正的状态。将测量状态值与加速度和速度的估计状态值进行比较可以提供检测机制。

图形子图（211），plot（t，ysig（1，:)  -  ysigest（1，:)）标题（'速度测量错误'）子平板（212），plot（t，ysig（2，：） -  ysigest（2，:)）标题（'加速度测量误差'）

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

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

图，子图（211），plot（t，fkur（1，:)）标题（'速度错误kurtosis')次要情节(212)、图(t, fKur(2:))标题('加速误差kurtosis'）

图包含2个轴。具有标题速度误差kurtosis的轴1包含型号的物体。带有标题加速度误差Kurtosis的轴2包含类型线的物体。

忽略前4秒估计估计仍然会聚并收集数据，误差的Kurtosis是相对恒定的，并且在3左右的次要变化（高斯分布的预期峰值值）。因此，错误统计信息不能用于自动检测此应用程序中的摩擦变化。在本申请中使用误差的kurtosis也很困难，因为过滤器正在调整和不断将误差驱动到零，只提供一个短时间窗口，其中错误分布不同于零。

因此，在本申请中，使用估计摩擦的变化提供了在电机中自动检测故障的最佳方法。来自已知无故障数据的摩擦估计（均值和标准偏差）为摩擦提供了期望的界限，并且何时违反这些界限时易于检测。以下绘图突出显示该故障检测方法。

图绘图（t，xsigest（2，：），[t nan t]，[fmean + 3 * fstd，nan，fmean-3 * fstd]标题（'摩擦变化检测'）传奇（“估计摩擦”，“无过错摩擦边界”）轴（[0 20 -10 20]）网格上

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

总结

此示例显示了如何使用扩展的卡尔曼滤波器来估计简单的直流电机中的摩擦，并使用摩擦估计进行故障检测。

使用扩展卡尔曼滤波器进行故障检测

汽车模型

创建扩展的卡尔曼过滤器

模拟

扩展卡尔曼过滤性能

总结

相关主题

预测维护工具箱文档

金宝app

用Matlab预测维护