基于残差分析的离心泵故障诊断

这个示例使用:

打开生活的脚本

这个例子展示了一种基于模型奇偶方程的方法来检测和诊断泵系统中发生的不同类型的故障。该示例扩展了基于稳态实验的离心泵故障诊断适用于泵跨多个工况运行的情况。

该实例遵循Rolf Isermann[1]所著的《故障诊断应用》一书中的离心泵分析。它使用了系统识别工具箱™、统计和机器学习工具箱™、控制系统工具箱™和Simulink™的功能，不需要预测性维护工具箱™。金宝app

多转速泵运行-残差分析诊断

如果泵在快速变化或更大的速度范围内运行，稳态泵扬程和扭矩方程不能产生准确的结果。摩擦和其他损失可能变得很重要，模型的参数显示出对速度的依赖。在这种情况下，一个广泛适用的方法是创建行为的黑盒模型。这些模型的参数不需要具有物理意义。该模型被用作模拟已知行为的设备。从相应的测量信号中减去模型的输出进行计算残差．残差的性质，如它们的均值，方差和幂被用来区分正常和错误的操作。

利用静态泵头方程和动态泵管方程，可以计算出图中所示的4个残差。

模型包含以下组件:

静态泵模型: $Δ p_{}^{ˆ} （ t ）＝ θ_{1} ω^{2} （ t ） + θ_{2} ω （ t ）$
管道动态模型: $问_{}^{ˆ} （ t ）＝ θ_{3.} + θ_{4} \sqrt{Δ p （ t ）} + θ_{5} 问_{}^{ˆ} （ t - 1 ）$
动态泵管模型: ${问_{}^{ˆ}}_{}^{ˆ} （ t ）＝ θ_{3.} + θ_{4} \sqrt{{Δ p （ t ）}_{}^{ˆ}} + θ_{5} {问_{}^{ˆ}}_{}^{ˆ} （ t - 1 ）$
动态逆泵模型: ${米_{}^{ˆ}}_{米 o t o r} （ t ）＝ θ_{6} + θ_{7} ω （ t ） + θ_{8} ω （ t - 1 ） + θ_{9} ω^{2} （ t ） + θ_{10} {米_{}^{ˆ}}_{米 o t o r} （ t - 1 ）$

模型参数 $θ_{1} ，．．．， θ_{10}$ 显示对泵速的依赖。在这个例子中，计算参数的分段线性逼近。将操作区域划分为3个区域:

$ω \leq 900 R P 米$
$900 < ω \leq 1500 R P 米$
$ω > 1500 R P 米$

一个健康的泵在一个闭环控制器的闭环中运行参考转速范围为0 - 3000rpm。参考输入是一个改进的PRBS信号。在10hz采样频率下采集电机转矩、泵转矩、转速和压力的测量值。加载测量到的信号并绘制参考和实际泵速(这些是大数据集，约20MB，可从MathWorks支持文件站点下载)。金宝app

url =“//www.tatmou.com/金宝appsupportfiles/predmaint/fault-diagnosis-of-centrifugal-pumps-using-residual-analysis/DynamicOperationData.mat”；websave (“DynamicOperationData.mat”url);负载DynamicOperationData图(t, RefSpeed, t, w) xlabel(“时间(s)”) ylabel (“泵转速(RPM)”)传说(“参考”，“实际”）

根据泵的转速范围定义操作制度。

I1 = w < = 900;第一运行状态I2 = w>900 & w<=1500;%秒运行状态I3 = w > 1500;第三种操作方式

模式识别

A.静压泵型号识别

估计的参数 $θ_{1}$ 和 $θ_{2}$ 在静态泵的方程中利用实测的泵速值 $ω （ t ）$ 和压差 $Δ p （ t ）$ 输入-输出数据。参见helper函数staticPumpEst执行这个估计。

th1 = 0 (3,1);th2 = 0 (3,1);dp =南(大小(dp));估计压差%[th1(1)， th2(1)， dpest(I1)] = staticPumpEst(w, dp, I1);形式1的估计[th1(2)， th2(2)， dpest(I2)] = staticPumpEst(w, dp, I2);* * * * * * * * * * * * * * * * * *[th1(3)， th2(3)， dpest(I3)] = staticPumpEst(w, dp, I3);* * * * * * * * * * * * * * * * * * * *Plot (t, dp, t, dpest)比较测量和预测压差包含(“时间(s)”) ylabel (‘\δP’)传说(“测量”，“估计”，“位置”，“最佳”)标题(“静态泵模型验证”）

B.动态管道模型识别

估计的参数 $θ_{3.}$ ， $θ_{4}$ 和 $θ_{5}$ 在管道流量方程中 $问_{}^{ˆ} （ t ）＝ θ_{3.} + θ_{4} \sqrt{Δ p （ t ）} + θ_{5} 问_{}^{ˆ} （ t - 1 ）$ ，使用测量的流量值 $问（ t ）$ 和压差 $Δ p （ t ）$ 输入-输出数据。参见helper函数dynamicPipeEst执行这个估计。

th3 = 0 (3,1);th4 = 0 (3,1);th5 = 0 (3,1);[th3 (1) th4 (1) th5 (1)] = dynamicPipeEst (dp Q I1);%(3)(4)(5)对于状态1的估计[th3 (2), th4 (2), th5 (2)] = dynamicPipeEst (dp Q I2);% ta3, ta4, ta5对状态2的估计[th3 (3), th4 (3), th5 (3)] = dynamicPipeEst (dp Q I3);% ta3, ta4, ta5对状态3的估计

与静态泵模型不同的是，动态管道模型表现出对流量值的动态依赖。在Simulink中，利用控制系统工具箱中的“LPV System”模块建立分段线性模型，对不同转速下的模型进行仿真。金宝app参见Simuli金宝appnk模型LPV_pump_pipe还有辅助函数simulatePumpPipeModel来执行模拟。

检查控制系统工具箱的可用性ControlsToolboxAvailable = ~ isempty(版本(“控制”) & &许可证(“测试”，“Control_Toolbox”）;如果ControlsToolboxAvailable%模拟动态管道模型。以测量的压力值作为输入Ts = t - t (1) (2);开关= 1(大小(w));开关(I2) = 2;开关(I3) = 3;UseEstimatedP = 0;Qest_pipe = simulatePumpPipeModel (Ts、th3 th4, th5);情节(t, Q t Qest_pipe)比较实测流量和预测流量其他的加载预保存的分段线性Simulink模型的仿真结果金宝app负载DynamicOperationDataQest_pipeTs = t - t (1) (2);情节(t, Q t Qest_pipe)结束包含(“时间(s)”) ylabel ('流量(Q)， m^3/s')传说(“测量”，“估计”，“位置”，“最佳”)标题(“动态管道模型验证”）

C.动态泵管模型识别

动态泵管模型使用与上述相同的参数( $θ_{3.} ， θ_{4} ， θ_{5}$ )，只是模型模拟需要使用估计的压差，而不是测量的压差。因此不需要新的身份证明。检查估算值 $θ_{3.} ， θ_{4} ， θ_{5}$ 给出了泵管动力学的良好再现。

如果ControlsToolboxAvailable UseEstimatedP = 1;Qest_pump_pipe = simulatePumpPipeModel (Ts、th3 th4, th5);情节(t, Q t Qest_pump_pipe)比较实测流量和预测流量其他的负载DynamicOperationDataQest_pump_pipe情节(t, Q t Qest_pump_pipe)结束包含(“时间(s)”) ylabel ('流量Q (m^3/s)')传说(“测量”，“估计”，“位置”，“最佳”)标题(“动态泵-管道模型验证”）

配合几乎与使用测量的压力值所获得的配合完全相同。

D.动态逆泵模型辨识

的参数 $θ_{6} ，．．．， θ_{10}$ 可以用类似的方式进行识别，通过将测量到的扭矩值回归到之前的扭矩和速度测量值。然而，一个完整的多速度仿真所得到的分段线性模型并不能很好地拟合数据。因此，尝试了一种不同的黑盒建模方法，包括识别一个非线性ARX模型与理性回归器来拟合数据。

%在550个样品中使用前300个样品进行鉴定N = 350;sys3 = identifyNonlinearARXModel (Mmot, w, Q, Ts, N)

sys3 =具有1个输出和2个输入的非线性ARX模型。变量y1 u1 u2 2的线性回归。自定义回归器:u1(t-2)^2自定义回归器:u1(t)*u2(t-2)自定义回归器:u2(t)^2所有回归器列表输出函数:线性函数样本时间:0.1秒状态:使用NLARX在时域数据上估计。拟合估计数据:49.2%(模拟焦点)FPE: 1798, MSE: 3392

Mmot_est = sim(sys3，[w Q]); / /用户名Mmot情节(t, t, Mmot_est)比较测量的和预测的电机转矩包含(“时间(s)”) ylabel (“电机转矩(Nm)”)传说(“测量”，“估计”，“位置”，“最佳”)标题(“反泵模型验证”）

残留一代

定义残留被测信号与相应模型产生的输出之间的差值。计算对应于四个模型组件的四个残差。

R1 = dp - dpest;r2 = Q - Qest_pipe;r3 = Q - Qest_pump_pipe;

在计算泵模型逆残差时，由于原残差方差较大，采用移动平均滤波器对模型输出进行平滑处理。

r4 = Mmot - movmean(Mmot_est，[1 5]);

对培训残留物的看法:

图subplot(221) plot(t,r1) ylabel('静压泵- r1') subplot(222) plot(t,r2)'动态管道- r2') subplot(223) plot(t,r3) ylabel('动力泵管道- r3')包含(“时间(s)”)子图(224)plot(t,r4) ylabel('动态逆泵- r4')包含(“时间(s)”）

残留物特征提取

残差是信号，从中提取合适的特征用于故障隔离。由于没有可用的参数信息，考虑纯粹由信号特性推导出来的特征，如信号的最大振幅或方差。

考虑使用PRBS输入实现对泵管系统进行的一组20个实验。以下每种模式均重复实验集:

健康的泵
故障1:间隙磨损
故障2:叶轮出口有少量沉积物
故障3:叶轮进口处有沉积物
故障4:叶轮出口磨料磨损
故障5:刀片断裂
错误6:空化
故障7:速度传感器偏置
故障8:流量计偏置
故障9:压力传感器偏置

加载实验数据

url =“//www.tatmou.com/金宝appsupportfiles/predmaint/fault-diagnosis-of-centrifugal-pumps-using-residual-analysis/MultiSpeedOperationData.mat”；websave (“MultiSpeedOperationData.mat”url);负载MultiSpeedOperationData%生成操作模式标签标签= {“健康”，“ClearanceGapWear”，“ImpellerOutletDeposit”，…“ImpellerInletDeposit”，“AbrasiveWear”，“BrokenBlade”，“空化”，“SpeedSensorBias”，…“FlowmeterBias”，“PressureSensorBias”}；

计算每个集合和每个操作模式的留数。这需要几分钟。因此，剩余数据被保存在数据文件中。运行helperComputeEnsembleResidues生成残差，如:

% HealthyR = helperComputeEnsembleResidues(HealthyEnsemble,Ts,sys3,th1,th2,th3,th4,th5);%正常数据残值

%加载预先保存的数据从"helperComputeEnsembleResidues"运行url =“//www.tatmou.com/金宝appsupportfiles/predmaint/fault-diagnosis-of-centrifugal-pumps-using-residual-analysis/Residuals.mat”；websave (“Residuals.mat”url);负载残差

残差的模分辨能力最大的特征是未知的。因此生成几个候选特征:均值，最大振幅，方差，峰度和每个残差的1范数。所有特性都使用健康集成中的值范围进行缩放。

CandidateFeatures = {@mean, @ (x)马克斯(abs (x)), @kurtosis, @var, @ (x)和(abs (x))};FeatureNames = {“的意思是”，“马克斯”，“峰度”，“方差”，“OneNorm”}；%从收集的每个故障模式残差生成特征表[HealthyFeature, MinMax] = helperGenerateFeatureTable(HealthyR, CandidateFeatures, FeatureNames);Fault1Feature = helperGenerateFeatureTable(Fault1R, CandidateFeatures, FeatureNames, MinMax);Fault2Feature = helperGenerateFeatureTable(Fault2R, CandidateFeatures, FeatureNames, MinMax);Fault3Feature = helperGenerateFeatureTable(Fault3R, CandidateFeatures, FeatureNames, MinMax);Fault4Feature = helperGenerateFeatureTable(Fault4R, CandidateFeatures, FeatureNames, MinMax);Fault5Feature = helperGenerateFeatureTable(Fault5R, CandidateFeatures, FeatureNames, MinMax);Fault6Feature = helperGenerateFeatureTable(Fault6R, CandidateFeatures, FeatureNames, MinMax);Fault7Feature = helperGenerateFeatureTable(Fault7R, CandidateFeatures, FeatureNames, MinMax);Fault8Feature = helperGenerateFeatureTable(Fault8R, CandidateFeatures, FeatureNames, MinMax);Fault9Feature = helperGenerateFeatureTable(Fault9R, CandidateFeatures, FeatureNames, MinMax);

每个特征表有20个特征(每个残差信号有5个特征)。每个表包含50个观察值(行)，每个实验一个。

N = 50;%在每个模式的实验次数FeatureTable = […[HealthyFeature (1: N,:), repmat(标签(1)[N, 1]));…[Fault1Feature (1: N,:), repmat(标签(2),[N, 1]));…[Fault2Feature (1: N,:), repmat(标签(3),[N, 1]));…[Fault3Feature (1: N,:), repmat(标签(4),[N, 1]));…[Fault4Feature (1: N,:), repmat(标签(5),[N, 1]));…[Fault5Feature (1: N,:), repmat(标签(6),[N, 1]));…[Fault6Feature (1: N,:), repmat(标签(7),[N, 1]));…[Fault7Feature (1: N,:), repmat(标签(8),[N, 1]));…[Fault8Feature (1: N,:), repmat(标签(9),[N, 1]));…[Fault9Feature (1: N,:), repmat(标签(10),(N, 1]]];FeatureTable.Properties。VariableNames{结束}=“条件”；%预览一些训练数据的样本disp(FeatureTable([2 13 37 49 61 62 73 85 102 120]，:))

Mean1非常刻薄的Mean3 Mean4 Max1 Max2 Max3 Max4 Kurtosis1 Kurtosis2 Kurtosis3 Kurtosis4 Variance1 Variance2 Variance3 Variance4 OneNorm1 OneNorm2 OneNorm3 OneNorm4条件  ________ _______ _______ _______ _______ ________ ________ _________ _________ _________ _________ __________ _________ _________ _________ _________ ________ ________ ________________ _________________________ 0.32049 0.6358 0.6358 0.74287 0.39187 0.53219 0.53219 0.041795 0.18957 0.089492 0.089489 0 0.45647 0.6263 0.6263 0.15642 0.56001 0.63541 0.63541 0.24258 {'Healthy'} 0.21975 0.19422 0.19422 0.83704 0 0.12761 0.12761 0.27644 0.18243 0.19644 0.19643 0.20856 0.033063 0.16772 0.16773 0.025119 0.057182 0.19344 0.19344 0.193440.063167 {' health '} 1 1 1 0.78284 0.94535 0.94535 0.9287 0.41371 0.45927 0.45926 0.5934 0.92332 1 1 0.80689 0.99783 1 1 0.9574 {' health '} -2.6545 0.90555 0.905550.86324 1.3037 0.7492 0.7492 0.97823 -0.055035 0.57019 0.57018 1.4412 1.4411 0.73128 0.73128 0.80573 3.2084 0.90542 0.90542 0.49257 {'ClearanceGapWear'} -0.89435 0.43384 0.43385 1.0744 0.82197 0.30254 0.30254 -0.022325 0.21411 0.098504{'ClearanceGapWear'} -1.0949 0.44616 0.44616 1.143 0.19719 0.19719 0.19719 -0.012055 -0.10819 0.32603 0.32603 -0.0033592 0.43341 0.12857 0.12857 0.038392 1.2238 0.44574 0.44574 0.093243{'ImpellerOutletDeposit'} -3.0312 0.9786 0.9786 0.75241 1.473 0.97839 0.97839 1.0343 0.0050647 0.4917 0.4917 0.89759 1.7529 0.9568 0.9568 0.8869 0.6536 0.978590.97859 - 0.62706 {' ImpellerOutletDeposit '}

分类器设计

A.使用散点图可视化模式可分离性

通过对特征的目视检查开始分析。为此，考虑故障1:间隙磨损。要查看哪些特征最适合检测此故障，请生成带有“Healthy”和“ClearanceGapWear”标签的特征散点图。

T = FeatureTable (:, 1:20);P = T.Variables;R = FeatureTable.Condition;I = strcmp (R,“健康”) | strcmp (R,“ClearanceGapWear”）;f =图;gplotmatrix (P(我:),[],R (I)) f.Position (3:4) = f.Position (3:4) * 1.5;

虽然不清晰可见，但第1列和第17列中的特性提供了最多的分隔。更仔细地分析这些特性。

f =图;名称= FeatureTable.Properties.VariableNames;J = [1 17];流('间隙故障的选定特征:%s\n'strjoin(名字(J),”、“））

间隙故障选择特征:Mean1, OneNorm1

gplotmatrix (P (I, 17 [1]), [], R (I))

从图中可以清楚地看出，特征Mean1和OneNorm1可用于分离健康模式和间隙间隙故障模式。可以对每种故障模式执行类似的分析。在所有情况下，都有可能找到一组特征来区分故障模式。因此检测错误的行为总是有可能的。然而,错隔离由于相同的特征会受到多种故障类型的影响，因此比较困难。例如，特征Mean1 (r1的平均值)和OneNorm1 (r1的1-norm)显示了许多故障类型的变化。有些故障(如传感器偏差)在许多特征上是可分离的，因此更容易分离。

对于三个传感器偏置故障，从散点图的人工检查中选取特征。

图;I = strcmp (R,“健康”) | strcmp (R,“PressureSensorBias”) | strcmp (R,“SpeedSensorBias”) | strcmp (R,“FlowmeterBias”）;J = [1 4 6 16 20];选定特征指数%流('传感器的选定特征'偏差:%s\n'strjoin(名字(J),”、“））

传感器偏差的选择特征:Mean1, Mean4, Max2, Variance4, OneNorm4

gplotmatrix (P (I, J), [], R (I))

所选特征的散点图表明，在一对或多对特征上可以区分4种模式。使用简化的特征集训练一个20个成员的Tree Bagger分类器，用于简化故障集(仅传感器偏差)。

rng默认的%的再现性Mdl = TreeBagger(20, FeatureTable(I，[J 21])，“条件”，…“OOBPrediction”，“上”，…“OOBPredictorImportance”，“上”）;图绘制(oobError (Mdl))包含(树木的数量) ylabel (“误分类概率”）

误分类误差小于3%。因此，有可能选择并使用较小的特征集来对特定的故障子集进行分类。

B.使用分类学习者App进行多类分类

上一节的重点是手工检查散点图，以减少特定故障类型的特征集。这种方法可能会很乏味，而且可能不会涵盖所有故障类型。您能否设计一个分类器，以更自动化的方式处理所有故障模式?在统计学和机器学习工具箱中有许多分类器。一个快速尝试和比较它们的方法是使用分类学习者应用程序。

启动Classification Learner App，并从工作区中选择FeatureTable作为新会话的工作数据。留出20%的数据(每个模式10个样本)用于抵抗验证。
选择所有下模型类型部分的主选项卡。然后按火车按钮。
在很短的时间内，大约训练了20个分类器。他们的准确性显示在他们的名字旁边的历史面板。线性支持向量机分类器表现最好，对保留样本的准确率达到86%。这个分类器在识别“ClearanceGapWear”时有一些困难，它有40%的时间将其分类为“ImpellerOutletDeposit”。
要获得性能的图形视图，请从主选项卡的plot部分打开一个Confusion Matrix图。该图显示了所选分类器(这里是线性支持向量机分类器)的性能。

将性能最好的分类器导出到工作空间，并使用它预测新的度量。

总结

设计良好的故障诊断策略可以减少服务停机时间和部件更换成本，从而节省运营成本。该策略得益于对操作机器动力学的良好知识，结合传感器测量来检测和隔离不同类型的故障。通过实例介绍了一种基于残差的离心泵故障诊断方法。当建模任务复杂且模型参数依赖于运行条件时，该方法是一种很好的替代基于参数估计和跟踪的方法。

基于残差的故障诊断方法包括以下步骤:

使用物理考虑或黑盒系统识别技术对系统的可测量输入和输出之间的动态进行建模。
计算残差作为测量和模型产生的信号之间的差异。残差可能需要进一步过滤以提高故障隔离性。
从每个残差信号中提取峰值幅值、功率、峰度等特征。
使用异常检测和分类技术对故障进行检测和分类。
并不是所有的残差和派生特征对每个故障都是敏感的。特征直方图和散点图的观点可以揭示哪些特征适合检测某一类型的故障。挑选特性并评估它们的性能以进行故障隔离的过程可能是一个迭代过程。

参考文献

Isermann,罗尔夫,故障诊断的应用。基于模型的状态监测:执行机构，驱动器，机械，工厂，传感器和容错系统，第1版，施普林格-弗拉格柏林海德堡，2011。

金宝app支持功能

静压泵方程参数估计

函数[x1, x2, dpest] = staticPumpEst(w, dp, I)%staticPumpEst在变速设置时估计的静态泵参数% I:所选操作区域的样本指数。w1 = [0;w(我)];dp1 = [0;dp (I)];(w1 R1 =。^ 2 w1);x = pinv (R1) * dp1;x1 = x (1);x2 = x (2);dp = R1(2:最终,:)* x;结束

动态管道参数估计

函数[x3, x4, x5, Qest] = dynamicPipeEst(dp, Q, I)%dynamicPipeEst动态管道参数估计在变化的速度设置% I:所选操作区域的样本指数。Q = Q(我);dp = dp (I);R1 = [0;问(1:end-1)];R2 = dp;R2 (R2 < 0) = 0;R2 = sqrt (R2);R = [ones(size(R2))， R2, R1];去除晶型外的样品2 =找到(我);j =找到(diff (ii) ~ = 1);R = R (2:,:);: R (j) = [];y = Q(2:结束);y (j) = [];x = R \ y;x3 = x (1);x4 = x (2);x5 = x (3); Qest = R*x;结束

利用LPV系统块对泵管模型进行多工况动态仿真。

函数q = simulatePumpPipeModel (Ts、th3 th4, th5)动态管道系统分段线性建模。% Ts:采样时间% w:泵转速% th1, th2, th3是估算的3种状态下的模型参数。这个功能需要控制系统工具箱。(魔法石,第1章=党卫军th5 (1) th4 (1) th5 (1) th4 (1) Ts);ss2 = ss (th5 (2), th4 (2), th5 (2), th4 (2) Ts);ss3 = ss (th5 (3), th4 (3), th5 (3), th4 (3), Ts);抵消=排列([th3 (1) th3 (2), th3(3)]的[3 2 1]);OP =结构(“地区”, (1 2 3) ');sys =猫(ss3 ss2 3魔法石,第1章);sys。SamplingGrid =运算;assignin (“基地”，“sys”sys) assignin (“基地”，“抵消”,抵消)mdl =“LPV_pump_pipe”；sim (mdl);q = logsout.get (“问”）;q = Qest.Values;q = Qest.Data;结束

建立了泵的逆动力学模型。

函数syse = identifyNonlinearARXModel (Mmot, w, Q, Ts, N)为2输入(w, Q)， 1输出(Mmot)数据识别一个非线性ARX模型。%的输入:% w:转速% Q:流量% Mmot:电机转矩% N:要使用的数据样本数量%输出:% syse:已识别的型号％该函数使用了系统识别工具箱中的NLARX估计器。Sys = idnlarx([2 2 1 0 1]， 1);＇＇，“CustomRegressors”, {“u1(2) ^ 2”，“u1 (t) * u2 (2)，“u2 (t) ^ 2”}）;数据= iddata(Mmot，[w Q]，Ts);选择= nlarxOptions;opt.Focus =“模拟”；opt.SearchOptions.MaxIterations = 500;syse = nlarx(数据(1:N)、系统选择);结束

基于残差分析的离心泵故障诊断

多转速泵运行-残差分析诊断

模式识别

残留一代

残留物特征提取

分类器设计

总结

参考文献

金宝app支持功能

相关的话题

预测性维护工具箱文档

金宝app

基于MATLAB的预测维护