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

打开生活的脚本

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

电机模型

电机被建模为一个带有阻尼系数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;％惯性t = 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) % %% u -电机扭矩输入% J -电机惯性% Ts -采样时间% %输出:% y -电机测量用元件[角速度;版权所有2016-2017 The MathWorks, Inc. % Extract data from input u = varargin{1};% Input 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} + （你_{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[q]=E[r]=0和协方差q = e [qq']和R = E (rr)．摩擦状态具有较高的过程噪声干扰。这反映了我们期望摩擦在电机正常运行时变化，并希望过滤器跟踪这种变化的事实。加速度和速度状态噪声很低，但速度和加速度测量相对有噪声。

指定过程噪声协方差。

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

指定测量噪声协方差。

r = [．..1e-40；．..%速度测量0 1的军医];%加速度测量

创建扩展卡尔曼滤波器

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

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

EKF = ExtendedKalmanFilter（．..@pdmMotorModelStateFcn,．..@PDM模型测量Cn，．..x0，．..“StateCovariance”, (1 0;0 1000),．..[1 0 0;0 1 0;00 100]，…'processnoise'问,．..“MeasurementNoise”R．..'stateTransitionJacobianfcn'，@pdmMotorModelStateJacobianFcn，．..'measurementjacobianfcn'，@pdmMotorModelMeasJacobianFcn）；

扩展卡尔曼滤波器将先前定义的状态转移和测量函数作为输入参数。初始状态值x0，初始状态协方差，过程和测量噪声协方差也是扩展卡尔曼滤波器的输入。在这个例子中，精确的雅可比函数可以从状态转移函数推导出来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=Pdmotormodelmeasjacobianfcn（x，varargin）%Pdmotormodelmeasjacobianfcn%%电机模型测量方程的雅可比矩阵。有关模型方程，请参见%Pdmotormodelmeasurementfcn。%%Jac=Pdmotormodelmeasjacobianfcn（x，u，J，Ts）%%输入：%x-元件状态[角速度；摩擦]%u-电机扭矩输入%J-电机惯性%Ts-采样时间%%输出：%Jac-测量雅可比矩阵在x%%版权所有2016-2017 MathWorks，Inc.%系统参数J=varargin{2}；%System innertia%Ja可比矩阵J=[…10；-x（2）/J-x（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）;估计摩擦力的峰度百分比fChanged = false (nt);%表示摩擦变化检测的标志

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

rng(“默认”）;Qv =胆固醇(Q);工艺噪声的标准偏差Qv（end）=1e-2；摩擦噪声小房车=胆固醇(R);测量噪声的％标准偏差

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

为了ct = 1:元素个数(t)型号输出更新y = pdmMotorModelMeasurementFcn (x0, u (ct), J, Ts);y = y +房车* 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=最大值（1，ct-400）：最大值（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）；%使用预期的摩擦平均值和标准偏差来检测%摩擦的变化。est摩擦=平均值（xSigEst（2，最大值（1，ct-10）：ct））；fChanged（ct）=（est摩擦>fMean（ct）+3*fSTD（ct））| |（est摩擦结尾如果fChanged (ct) & & ~ fChanged (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））;．..idx峰度(ySigEst (2) -ySig (2, idx)));结尾

扩展卡尔曼滤波性能

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

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

图中包含2个轴对象。标题为Motor Outputs的轴对象1包含2个类型为line的对象。这些物体代表测量的角速度，测量的角加速度。标题为“电机输入-扭矩”的轴对象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个类型为line的对象。这些对象分别表示真值、估计值、置信区间。轴对象2，标题为电机状态-摩擦包含3个类型为线的对象。

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

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

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

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

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

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

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

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

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

忽略估计器仍在收敛且正在收集数据时的前4秒，误差的峰度相对恒定，约为3（高斯分布的预期峰度值）。因此，在此应用中，误差统计信息不能用于自动检测摩擦变化。在该应用中，使用误差峰度也很困难，因为滤波器正在适应并不断地将误差驱动到零，仅给出一个误差分布不同于零的短时间窗口。

因此，在本应用中，使用估计摩擦力的变化提供了自动检测电机故障的最佳方法。摩擦力估计值（平均值和标准偏差）根据已知的无故障数据，可以提供摩擦的预期界限，并且很容易在违反这些界限时进行检测。下图突出显示了这种故障检测方法。

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

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

总结

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

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

电机模型

创建扩展卡尔曼滤波器

模拟

扩展卡尔曼滤波性能

总结

相关的话题

预测性维护工具箱文档

金宝app

基于MATLAB的预测维修