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

这个例子使用了:

打开实时脚本

本例展示了一种基于模型奇偶方程的方法，用于检测和诊断泵系统中发生的不同类型的故障。中的技术进行了扩展基于稳态实验的离心泵故障诊断到泵运行跨越多种工况的情况。

该示例遵循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)“时间(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);dest = nan(大小(dp));估计压差%[th1(1)， th2(1)， dest (I1)] = staticPumpEst(w, dp, I1);区域1的估计值[th1(2)， th2(2)， dest (I2)] = staticPumpEst(w, dp, I2);区域2的估计值[th1(3)， th2(3)， dest (I3)] = staticPumpEst(w, dp, I3);区域3的估计值Plot (t, dp, t, dest)%比较实测和预测的压差包含(“时间(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);% 3 4 5在区域2中的估计值[th3 (3), th4 (3), th5 (3)] = dynamicPipeEst (dp Q I3);% 3, 4, 5对区域3的估计

与静态泵模型不同，动态管道模型显示了流量值的动态依赖性。为了模拟不同转速下的模型，在Simulink中使用Control System Toolbox中的“LPV System”块创建了分段线性模型。金宝app参见Simuli金宝appnk模型LPV_pump_pipe还有辅助函数simulatePumpPipeModel执行模拟。

检查控制系统工具箱是否可用ControlsToolboxAvailable = ~isempty(ver(“控制”) &&许可证(“测试”，“Control_Toolbox”）;如果ControlsToolboxAvailable模拟动态管道模型。使用测量的压力值作为输入Ts = t(2)-t(1);Switch = ones(size(w));开关(I2) = 2;开关(I3) = 3;UseEstimatedP = 0;Qest_pipe = simulatePumpPipeModel(Ts,th3,th4,th5);情节(t, Q t Qest_pipe)比较实测流量和预测流量其他的负载预保存的分段线性Simulink模型仿真结果金宝app负载DynamicOperationDataQest_pipeTs = t(2)-t(1);情节(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模型输入:u1, u2输出:y1标准回归量对应阶数:na = [2] nb = [2 1] nk =[0 1]自定义回归量:u1(t-2)^2 u1(t)*u2(t-2) u2(t)^2非线性回归量:无模型输出在回归量中是线性的。采样时间:0.1秒状态:使用时域数据的NLARX估计。拟合估计数据:50.55%(模拟焦点)FPE: 1759, MSE: 3214

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

残留一代

定义残留一个模型的测量信号与相应的模型输出之间的差值。计算四个模型分量对应的四个残差。

R1 = dp - dest;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)动态泵管- r3)包含(“时间(s)”) subplot(224) subplot(t,r4)动态逆泵- r4)包含(“时间(s)”）

残留特征提取

残差是信号，从中提取合适的特征用于故障隔离。由于没有参数信息可用，所以考虑纯粹从信号属性(如信号的最大幅度或方差)中派生的特征。

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

健康的泵
故障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 = helpcomputeensembler余量(HealthyEnsemble,Ts,sys3,th1,th2,th3,th4,th5);%健康数据残差

从“helpercomputeensemblerremainder”运行加载预先保存的数据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))};featurename = {“的意思是”，“马克斯”，“峰度”，“方差”，“OneNorm”}；%由收集的各故障模式残差生成特征表[HealthyFeature, MinMax] = helperGenerateFeatureTable(HealthyR, CandidateFeatures, featurename);Fault1Feature = helperGenerateFeatureTable(Fault1R, CandidateFeatures, featurename, MinMax);Fault2Feature = helperGenerateFeatureTable(Fault2R, CandidateFeatures, featurename, MinMax);Fault3Feature = helperGenerateFeatureTable(Fault3R, CandidateFeatures, featurename, MinMax);Fault4Feature = helperGenerateFeatureTable(Fault4R, CandidateFeatures, featurename, MinMax);Fault5Feature = helperGenerateFeatureTable(Fault5R, CandidateFeatures, featurename, MinMax);Fault6Feature = helperGenerateFeatureTable(Fault6R, CandidateFeatures, featurename, MinMax);Fault7Feature = helperGenerateFeatureTable(Fault7R, CandidateFeatures, featurename, MinMax);Fault8Feature = helperGenerateFeatureTable(Fault8R, CandidateFeatures, featurename, MinMax);Fault9Feature = helperGenerateFeatureTable(Fault9R, CandidateFeatures, featurename, 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.45647 0.6263 0.6263 0.53219 0.041795 0.18957 0.089492 0.089489 00.15642 0.56001 0.63541 0.63541 0.24258{'健康'}0.21975 0.19422 0.19422 0.83704 0 0.12761 0.12761 0.27644 0.19643 0.20856 0.033063 0.16761 0.16773 0.025119 0.19344 0.063167{'健康'}0.1847 0.25136 0.13929 0.1393 0.084162 0.72346 0.091488 0.091485 0.052552 0.063757 0.18371 0.18371 0.023845 0.10671 0.25065 0.25065 0.04848{'健康'}11 10 0.78284 0.94535 0.94535 0.75934 0.92332 11 0.80689 0.9978311 0.9574 {'Healthy'} -2.6545 0.90555 0.90555 0.86324 1.3037 0.7492 0.7492 0.97823 -0.055035 0.57019 0.57018 1.4412 1.4411 0.73128 0.90542 0.80573 3.2084 0.90585 1.0744 0.82197 0.30254 0.30254 -0.022325 0.21411 0.098504 0.0985 -0.010587 0.55959 0.29554 0.29554 0.066693 1.0432 0.3326 0.43326 0.13785 {'ClearanceGapWear'} -1.2149 0.53579 0.53579 1.0899 0.4739 0.47339 0.47339 0.29539 0.048137 0.17864 0.796480.48658 0.48658 0.25686 1.4066 0.5353 0.5353 0.26663 {'ClearanceGapWear'} -1.0949 0.44616 0.44616 1.143 0.56693 0.19719 0.19719 - 0.032603 0.32603 -0.0033592 0.43341 0.12857 0.12857 0.038392 0.2238 0.44574 0.093243 {'ClearanceGapWear'} -0.1844 0.16651 0.16651 0.87747 0.25597 0.080265 0.080265 0.49715 0.5019 0.29939 0.29938 0.5338 0.072397 0.058024 0.058025 0.1453 0.20263 0.16566 0.16566 0.14216 {'ImpellerOutletDeposit'} -3.0312 0.9786 0.75241 1.473 0.978390.97839 1.0343 0.0050647 0.4917 0.4917 0.89759 1.7529 0.9568 0.9568 0.8869 3.6536 0.97859 0.97859 0.62706 {'ImpellerOutletDeposit'}

分类器设计

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

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

T = FeatureTable(:，1:20);P = t变量;R = featutable . 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 =数字;名称= featutable . properties . variablenames;J = [1 17];流(间隙故障选择特征:%s\nstrjoin(名字(J),”、“）)

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

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

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

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

图;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，[21])，“条件”，.．.“OOBPrediction”，“上”，.．.“OOBPredictorImportance”，“上”）;图(oobError(Mdl)) xlabel(“树的数量”) ylabel (“误分类概率”）

误分类误差小于3%。因此，可以选择和使用更小的特征集对故障的某个子集进行分类。

B.使用Classification Learner App进行多类分类

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

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

将性能最好的分类器导出到工作区，并将其用于对新度量的预测。

总结

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

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

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

参考文献

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

金宝app支持功能

静态泵方程参数估计

函数[x1, x2, dest] = staticPumpEst(w, dp, I)静态泵参数估计在变化的速度设置% I:所选操作区域的样本指数。W1 = [0;w(我)];Dp1 = [0;dp (I)];R1 = [w1.]^ 2 w1);x = pinv(R1)*dp1;X1 = x(1);X2 = x(2);dest = R1(2:end，:)*x;结束

动态管道参数估计

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

基于LPV系统块的泵管模型动态多工况仿真。

函数Qest = simulatePumpPipeModel(Ts,th3,th4,th5)动态管道系统分段线性建模。% Ts:采样时间% w:泵转速% th1, th2, th3为3种状态下的估计模型参数。此功能需要控制系统工具箱。ss1 = ss(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 = struct(“地区”，[1 2 3]');Sys = cat(3,ss1,ss2,ss3);sys。SamplingGrid = OP;assignin (“基地”，“sys”sys) assignin (“基地”，“抵消”，offset) MDL =“LPV_pump_pipe”；sim (mdl);Qest = logsout.get(“问”）;Qest = Qest. values;Qest = 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]，”，“CustomRegressors”, {“u1(2) ^ 2”，“u1 (t) * u2 (2)，“u2 (t) ^ 2”});data = iddata(Mmot，[w Q]，Ts);opt = nlarxOptions;opt.Focus =“模拟”；opt.SearchOptions.MaxIterations = 500;syse = nlarx(data(1:N)，sys,opt);结束

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

多速泵运行-残余分析诊断

模式识别

残留一代

残留特征提取

分类器设计

总结

参考文献

金宝app支持功能

相关的话题

预测性维护工具箱文档

金宝app

基于MATLAB的预测性维护