基于模拟数据的多类故障检测
这个示例展示了如何使用Simulink模型来生成故障和健康数据。万博1manbetx这些数据被用来开发一个多类分类器来检测不同的故障组合。该实例使用了三缸往复泵模型，包括泄漏、堵塞和轴承故障。
设置模型
这个示例使用了许多存储在zip文件中的支持文件万博1manbetx。解压缩文件以访问支持文件，加载模型参数，并创建往复泵库。万博1manbetx
如果~ (“+机械液压力”，“dir”)解压缩('pdmrecippump_万博1manbetxsupportingfiles.zip'）结束％加载参数pdmRecipPump_Parameters％泵CAT_Pump_1051_DataFile_imported%计算机辅助设计如果需要，％创建Simscape库如果存在('mech_hydro_forcesps_lib'，“文件”) ~ = 4 ssc_build机械液压泵结束
往复泵模型
往复泵由电动机、泵壳、泵曲柄和泵柱塞组成。
mdl ='pdmrecippump'; 开放式系统（mdl）
(mdl open_system (,“/泵”])
所述泵型号配置为模拟三种类型的故障;汽缸泄漏，入口堵塞，轴承摩擦增加。这些错误被参数化为工作空间变量，并通过泵块对话框进行配置。
模拟故障和正常数据
为这三种故障类型中的每一种创建一个表示故障严重程度的值数组，范围从无故障到严重故障。
%定义故障参数变化numparvalues = 10;Leak_area_set_factor = linspace（0.00,0.036，numparvalues）;leak_area_set = leak_area_set_factor * trp_par.check_valve.in.max_area;LEAD_AREA_ET = MAX（LEAD_AREA_SET，1E-9）;%泄漏区域不能为0blockinfactor_set=linspace（0.8,0.53，numParValues）；承载因子_集=linspace（0,6e-4，数值）；
泵模型被配置为包括噪声，从而运行具有相同故障参数值的模型，将导致不同的仿真输出。这对于开发分类器是有用的，因为它意味着可以有多个模拟结果对于相同的故障状况和严重性。要为此类结果配置模拟，请创建故障参数值的向量，其中值表示没有故障，单个故障，两个故障的组合以及三个故障的组合。对于每个组（无故障，单个故障等），从上面定义的故障参数值创建125个故障值组合。这为此提供了1000个故障参数值组合。
nPerGroup = 125;%每个故障组中的元件数量无故障模拟leakArea = repmat (leak_area_set (1) nPerGroup, 1);blockingFactor = repmat (blockinfactor_set (1) nPerGroup, 1);bearingFactor = repmat (bearingfactor_set (1) nPerGroup, 1);%单故障模拟idx =装天花板(10 *兰德(nPerGroup, 1));leakArea = [leakArea;leak_area_set (idx)];blockingFactor = [blockingFactor; repmat (blockinfactor_set (1) nPerGroup, 1)];bearingFactor = [bearingFactor; repmat (bearingfactor_set (1) nPerGroup, 1)];idx =装天花板(10 *兰德(nPerGroup, 1));leakArea = [leakArea;nPerGroup repmat (leak_area_set (1), 1)];blockingFactor = [blockingFactor; blockinfactor_set (idx)];bearingFactor = [bearingFactor; repmat (bearingfactor_set (1) nPerGroup, 1)]; idx = ceil(10*rand(nPerGroup,1)); leakArea = [leakArea; repmat(leak_area_set(1),nPerGroup,1)]; blockingFactor = [blockingFactor;repmat(blockinfactor_set(1),nPerGroup,1)]; bearingFactor = [bearingFactor;bearingfactor_set(idx)'];双故障模拟idxA =装天花板(10 *兰德(nPerGroup, 1));idxB =装天花板(10 *兰德(nPerGroup, 1));leakArea = [leakArea;leak_area_set (idxA)];blockingFactor = [blockingFactor; blockinfactor_set (idxB)];bearingFactor = [bearingFactor; repmat (bearingfactor_set (1) nPerGroup, 1)];idxA =装天花板(10 *兰德(nPerGroup, 1));idxB =装天花板(10 *兰德(nPerGroup, 1));leakArea = [leakArea;leak_area_set (idxA)]; blockingFactor = [blockingFactor;repmat(blockinfactor_set(1),nPerGroup,1)]; bearingFactor = [bearingFactor;bearingfactor_set(idxB)']; idxA = ceil(10*rand(nPerGroup,1)); idxB = ceil(10*rand(nPerGroup,1)); leakArea = [leakArea; repmat(leak_area_set(1),nPerGroup,1)]; blockingFactor = [blockingFactor;blockinfactor_set(idxA)']; bearingFactor = [bearingFactor;bearingfactor_set(idxB)'];%三重故障模拟idxA =装天花板(10 *兰德(nPerGroup, 1));idxB =装天花板(10 *兰德(nPerGroup, 1));idxC =装天花板(10 *兰德(nPerGroup, 1));leakArea = [leakArea;leak_area_set (idxA)];blockingFactor = [blockingFactor; blockinfactor_set (idxB)];bearingFactor = [bearingFactor; bearingfactor_set (idxC)];
使用故障参数组合来创建万博1manbetx仿真软件。SimulationInput物体。对于每个模拟输入，确保设置不同的随机种子以生成不同的结果。
为simminput (ct) = Simulink.SimulationInput万博1manbetx(mdl);simInput (ct) = setVariable (simInput (ct),“泄漏\气缸\区域\ WKSP”渗漏面积（ct）；simInput（ct）=设置变量（simInput（ct），'block_in_factor_wksp'blockingFactor (ct));simInput (ct) = setVariable (simInput (ct),“bearing_fault_frict_WKSP”，承载因子（ct）；simInput（ct）=设置变量（simInput（ct），“noise_seed_offset_WKSP”，CT-1）;结束
使用generateSimulationEnsemble函数来运行由万博1manbetx仿真软件。SimulationInput上面定义的对象并将结果存储在本地子文件夹中。然后创建一个simulationEnsembleDatastore从存储的结果。
请注意，在标准桌面上并行运行这1000个模拟大约需要一个小时，并生成大约620MB的数据。为方便起见，提供了仅运行前10个模拟的选项。
％运行模拟并创建一个组合以管理仿真结果runAll = true;如果runall [确定，e] =生成（SimInput，FullFile（“。”，“数据”),“UseParallel”,真正的);其他的[ok，e]=生成模拟集成（模拟输入（1:10），完整文件(“。”，“数据”）);% #好< UNRCH >结束
[09-Apr-2018 09:01:38]正在检查并行池的可用性[2018年4月9日09:01:38]在并行工作程序上加载Simulink。。。分析并将文件传输给工人…完成[2018年4月9日09:万博1manbetx01:38]正在配置并行工作进程上的模拟缓存文件夹[2018年4月9日09:01:38]在平行工上运行SetupFcn[2018年4月9日09:01:39]平行作业工人的装载模型[2018年4月9日09:01:39]正在将模型中使用的基本工作空间变量传输给并行工作者[2018年4月9日09:01:41]运行模拟[2018年4月09日09:02:28]完成1000次模拟运行中的1次[2018年4月09:02:33]完成1000次模拟运行中的2次[2018年4月09:02:37]完成1000次模拟运行中的3次[2018年4月09:02:41]完成1000次模拟运行中的4次[2018年4月09:02:46]完成1000次模拟运行中的5次[2018年4月09:02:49]完成1000次模拟运行中的6次运行[09-Apr-2018 09:02:54]完成了1000次模拟运行中的7次[09-Apr-2018 09:02:58]完成了1000次模拟运行中的8次[09-Apr-2018 09:03:01]完成了1000次模拟运行中的9次。。。
ens=模拟安装数据存储（完整文件(“。”，“数据”）);
对仿真结果进行特征处理和提取
该模型配置为记录泵输出压力、输出流量、电机转速和电机电流。
数据变量
ans=8×1字符串数组"SimulationInput" "SimulationMetadata" "iMotor_meas" "pIn_meas" "pOut_meas" "qIn_meas" "qOut_meas" "wMotor_meas"
对于集合中的每个成员，对泵输出流量进行预处理，并根据泵输出流量计算工况指标。这些状态指标用于故障分类。为了进行预处理，要去掉输出流量的前0.8秒，因为这包含了模拟和泵启动的瞬态。作为预处理的一部分，计算输出流的功率谱，并使用SimulationInput返回故障变量的值。
配置集合，以便只读取兴趣的变量并调用预处理函数，该函数在本示例的最后定义。
ens.SelectedVariables = [“qOut_meas”，“SimulationInput”]; 重置（ens）data = read(ens)
data =1×2表qOut表示模拟投入量[2001×1时间表][1×1模拟投入量]万博1manbetx
(流、flowP flowF faultValues] =预处理(数据);
绘制流量和流量频谱。绘制的数据是为无故障的情况。
%标称数字次要情节(211);情节(flow.Time flow.Data);次要情节(212);semilogx (flowF pow2db (flowP));包含(“赫兹”）
流程谱显示在预期频率下的共振峰。具体而言，泵电机速度为950 RPM，或15.833 Hz，由于泵有3个气缸，预计流量将在3×15.833 Hz，47.5 Hz的基础上具有基础，以及47.5 Hz的倍数的谐波。流程谱清楚地显示了预期的共振峰。泵的一个气缸中的缺陷将导致泵电机速度的共振，15.833 Hz及其谐波。
流谱和慢信号给出了一些可能的条件指标。具体来说，是常见的信号统计量，如均值、方差等，以及频谱特征。计算与预期谐波有关的频谱状况指标，例如峰值幅度的频率、15.833 Hz左右的能量、47.5 Hz左右的能量、100 Hz以上的能量。并计算了谱峰的频率。
配置条件指标的数据变量和故障变量值的条件变量。然后调用extractCI函数来计算特征，并使用writeToLastMemberRead命令，将特性和故障变量值添加到集成。的extractCI函数在本例末尾定义。
ens.DataVariables = [ens.DataVariables;．..“fpeak”；“犁”；“pMid”；“pHigh”；“pKurtosis”；．..“qMean”；“qVar”；“qSkewness”；“qkurtosis”；．..“qPeak2Peak”；“qCrest”；“则”；“qMAD”；“qCSRange”]; ens.ConditionVariables=[“泄漏跳”，“BlockingFault”，“BearingFault”];壮举= extractCI(流、flowP flowF);datatwrite = [faultValues, feat];writeToLastMemberRead(实体,dataToWrite {:})
上述代码预处理并计算模拟集合第一个成员的条件指示器。使用集合对集合中的所有成员重复此操作hasdata指挥部。为了了解不同故障条件下的模拟结果，请绘制集合的每一百个元素。
%带有标称和故障的图形图，子图（211）;lflow = plot（flow.time，flow.data，“线宽”2);次要情节(212);lFlowP = semilogx (flowF pow2db (flowP),“线宽”,2); xlabel(“赫兹”)ct=1；lColors=get（lFlow.Parent，“ColorOrder”）;lIdx = 2;preprocess . %循环系统中所有成员%并计算每个成员的特征尽管hasdata(实体)读取成员数据data =阅读(ens);%对成员数据进行预处理和提取特征(流、flowP flowF faultValues] =预处理(数据);壮举= extractCI(流、flowP flowF);%将提取的特征值添加到成员数据datatwrite = [faultValues, feat];writeToLastMemberRead(实体,dataToWrite {:})%绘制每100个成员的成员信号和频谱如果Mod (ct,100) == 0 line('父母', lFlow。父母,'xdata'，flow.Time，“伊达塔”,流。数据,．..“颜色”lColors (lIdx:));线('父母', lFlowP。父母,'xdata'flowF,“伊达塔”pow2db (flowP),．..“颜色”lColors (lIdx:));如果lIdx==大小（l颜色，1）lIdx=1；其他的lIdx=lIdx+1；结束结束ct=ct+1；结束
请注意，在不同的故障条件和严重程度下，频谱包含在预期频率的谐波。
泵故障的检测和分类
前一节对仿真集合中所有成员的流量信号进行预处理和条件指标计算，对应不同故障组合和严重程度的仿真结果。状态指示器可用于从泵流量信号中检测和分类泵故障。
配置仿真集成以读取条件指示器，并使用tall和GARGET命令将所有条件指示器和故障变量值加载到内存中
%获取数据以设计分类器。重置(ens) ens. selectedvariables = [．..“fpeak”，“犁”，“pMid”，“pHigh”，“pKurtosis”，．..“qMean”，“qVar”，“qSkewness”，“qkurtosis”，．..“qPeak2Peak”，“qCrest”，“则”，“qMAD”，“qCSRange”，．..“泄漏跳”，“BlockingFault”，“BearingFault”];IdxlastFeature = 14;%加载状态指示器数据到内存中数据=聚集（高（ens））;
使用“local”配置文件启动并行池(parpool)…保留id为1 2的作业，因为它们包含崩溃转储文件。你可以使用'delete(myCluster.Jobs)'来删除所有用profile local创建的作业。要创建'myCluster'，请使用'myCluster = parcluster('local')'。连接到6个工人。使用Parallel Pool 'local'计算tall表达式
数据（1:10，:)
ans =表10×17fPeak犁pMid pHigh pKurtosis qMean qVar qSkewness qKurtosis qPeak2Peak查收我们则qMAD qCSRange LeakFault BlockingFault BearingFault ______ _______ ______ ______ _________ ______ ______ _________ _________ __________ ______ ______ ______ ________ _________ _____________ ____________ 43.909 0.86472 117.63 18.874 276.49 35.572 7.5242 -0.728322.7738 13.835 1.1494 35.677 2.2326 42690 0 1 e-09 0.8 43.909 0.44477 125.92 18.899 12.417 35.576 7.869 -0.7094 2.6338 13.335 1.1449 35.686 2.3204 42697 0 1 e-09 0.8 43.909 1.1782 137.99 17.526 11.589 35.573 7.4367 -0.72208 2.7136 12.641 1.1395 35.678 2.2407 42695 1 e-09 0.8 0 44.151 156.74 173.84 21.073 199.5 33.768 12.466 -0.30256 2.4782 17.4461.2138 33.952 2.8582 40518 2.4 e-06 0.8 0 43.848 0.71756 110.92 18.579 197.02 35.563 7.5781 -0.72377 2.793 14.14 1.1504 35.669 2.2671 42682 1 e-09 0.8 0 43.909 0.43673 119.56 20.003 11.589 35.57 7.5028 -0.74797 2.7913 13.833 1.1551 35.676 2.2442 42689 0 1 e-09 0.8 43.788 0.31617 135.3 19.724 476.82 35.568 7.4406 -0.70964 2.6884 14.685 1.1473 35.67342687 - 2.2392 1 e-09 0.8 0 43.848 0.72747 121.63 19.733 11.589 35.523 7.791 -0.72736 2.7864 14.043 1.1469 35.633 2.2722 42633 1 e-09 0.8 0 43.848 0.62777 128.85 19.244 11.589 35.541 7.5698 -0.6953 2.6942 13.451 1.1415 35.647 2.2603 42654 0 1 e-09 0.8 43.848 0.4631 134.83 18.918 12.417 35.561 7.8607 -0.68417 2.6664 13.885 1.1504 35.671 2.3078 426811 e-09 0.8 0
可以将每个集成成员(数据表中的行)的故障变量值转换为故障标志，并将故障标志组合为单个标志，以捕获每个成员的不同故障状态。
%将错误变量值转换为标志data.leakflag = data.leakfault> 1e-6;data.blockingflag = data.blockingfault <0.8;data.bearingflag = data.bearingfault> 0;data.combinedflag = data.leakflag + 2 * data.blockingflag + 4 * data.bearingflag;
创建一个分类器，将条件指示器作为输入，并返回组合故障标志。训练使用二阶多项式核的支持向量机。使用万博1manbetxcvpartition命令将集合成员划分为一个用于训练的集合和一个用于验证的集合。
rng(“默认”）%的再现性预测=数据(:,1:idxLastFeature);响应= data.CombinedFlag;本量利= cvpartition(大小(预测,1),“KFold”5);％创建和培训分类器template=templateSVM(．..'骨箱'，“多项式”，．..“PolynomialOrder”2,．..“KernelScale”，“汽车”，．..“BoxConstraint”,1．..“标准化”，对）；组合分类器=fitcecoc(．..预测(cvp.training (1):)．..响应(cvp.training (1):)．..“学习者”，模板，．..“编码”，'Onevsone'，．..“类名”, (0;1;2;3;4;5;6;7);
使用验证数据检查经过训练的分类器的性能，并将结果绘制在混淆图上。
通过计算和绘制混淆矩阵来检查性能actualValue =响应(cvp.test (1):);predictedValue = predict(combinedClassifier, predictors(cvp.test(1)，:));confdata = confusionmat (actualValue predictedValue);图，标签= {“没有”，“泄漏”，“阻止”，“泄露”&阻塞，“轴承”，．..'轴承和泄漏'，'轴承和阻挡'，“全部”};h =热图(confdata,．..'ylabel'，“实际泄漏故障”，．..“YDisplayLabels”、标签．..“XLabel”，'预测过错'，．..'xdisplaylabels'、标签．..“ColorbarVisible”，“关闭”）;
混淆图显示了每个故障组合正确预测故障组合的次数(图的对角线项)和错误预测故障组合的次数(非对角线项)。
混淆图显示分类器未正确分类某些故障条件（非对角项）。但是，无故障状态的预测是正确的。在一些地方，当存在故障（第一列）时，会预测无故障情况，否则会预测故障，尽管它可能不完全是正确的故障情况。总的来说，验证准确率为84%，预测故障的准确率为98%。
％计算分类器的整体准确性sum（diag（confdata））/sum（confdata（：）
ans=0.6150
%计算分类器在预测时的准确率％有错1-sum (confdata(2:最终,1))/笔(confdata (:))
ans = 0.9450
检查没有预测到故障但确实存在故障的情况。首先在验证数据中找到实际故障为阻塞故障但预测为无故障的情况。
vData =数据(cvp.test (1):);b1 = (actualValue==2) & (predictedValue==0);fData = vData (b1,十五17)
fData =11×3表泄漏故障阻塞故障轴承故障0.7701E-090.7701E-090.7701E-090.7701E-090.7701E-090.7701E-090.7701E-090.7701E-090.074
在验证数据中找出实际故障是轴承故障但预测没有故障的情况。
b2=（实际值==4）和（预测值==0）；vData（b2,15:17）
ans=0×3空表
检查没有故障预测的情况，但存在故障的情况显示，当堵塞故障值为0.77的阻塞故障靠近其标称值0.8时，或者轴承故障值接近其标称值0.用小阻塞故障值绘制频谱，并与无故障情况进行比较，显示光谱非常相似的使检测难以。重新培训分类器，但包括0.77的阻塞值，因为非故障条件会显着提高故障检测器的性能。或者，使用额外的泵测量可以提供更多信息并提高检测小阻挡故障的能力。
％配置集合以仅读取流量和故障变量值ens.SelectedVariables = [“qOut_meas”，“泄漏跳”，“BlockingFault”，“BearingFault”]; 重置（ens）%将集合成员数据加载到内存中数据=聚集（高（ens））;
使用并行池“local”计算tall表达式：-通过1/1：在38秒内完成计算在38秒内完成
%查找错误预测的成员，然后％计算其功率谱idx =．..data.BlockingFault==fData.BlockingFault（1）&．..数据。LeakFault == fData.LeakFault(1) &．..数据。BearingFault = = fData.BearingFault (1);flow1 =数据(idx, 1);flow1 = flow1.qOut_meas {1};[flow1P, flow1F] = pspectrum (flow1);%查找一个没有任何错误的成员idx =．..数据。BlockingFault == 0.8 &．..data.LeakFault==1e-9&．..数据。BearingFault = = 0;flow2 =数据(idx, 1);flow2 = flow2.qOut_meas {1};[flow2P, flow2F] = pspectrum (flow2);%绘制功率谱半对数(．..flow1F pow2db (flow1P),．..flow2F pow2db (flow2P));包含(“赫兹”)传说('小阻挡过错'，'没有错'）
结论
该实例展示了如何利用Simulink模型对往复泵故障进行建模，仿真万博1manbetx该模型在不同故障组合和严重程度下，从泵输出流量中提取工况指标，并用工况指标训练分类器检测泵故障。该示例检查了使用分类器进行故障检测的性能，并指出小阻塞故障与无故障情况非常相似，无法可靠地检测。
万博1manbetx支持功能
功能[流量、流量谱、流量频率、故障值]=预处理（数据）%辅助功能，用于预处理记录的往复泵数据。%移除流量信号的前0.8秒tmin =秒（0.8）;flow = data.qout_meas {1};flow = flow（flow.time> = tmin，:);flow.time = flow.time - flow.time（1）;%确保以统一的采样率进行流量采样流=调整时间(流,“普通”，“线性”，“步伐”秒(1 e - 3));%从流量中去除平均值，计算流量谱fA =流;足总。Data =。数据——意味着(fA.Data);[flowSpectrum, flowFrequencies] = pspectrum(足总,“FrequencyLimits”[250]);%从SimulationInput中找到故障变量的值思敏= data.SimulationInput {1};var = {simin.Variables.Name};idx = strcmp (var,“泄漏\气缸\区域\ WKSP”）;LeakFault = simin.Variables (idx) value;idx = strcmp (var,'block_in_factor_wksp'）;BlockingFault = simin.Variables (idx) value;idx = strcmp (var,“bearing_fault_frict_WKSP”); BearingFault=相似变量（idx）.值；%收集单元格数组中的故障值FaultValue={．..“LeakFault”，泄漏故障，．..“BlockingFault”，阻塞事实，．..'autingfault'，承载断层}；结束功能CI = ExtractCi（流量，流量，流量）％辅助功能从流量信号提取条件指示灯%和频谱。%找出功率谱中峰值幅值的频率。pMax = max (flowP);fPeak = flowF (flowP = = pMax);％计算低频范围10-20 Hz的功率。>= 10 & > <= 20;犁=总和(flowP(纤毛刷));%计算中频范围40hz ~ 60hz的功率。法兰=流量>=40，流量<=60；pMid=总和（流量（fRange））；%计算高频范围> 100hz的功率。>= 100;pHigh =总和(flowP(纤毛刷));%找到光谱峰度峰的频率[PKUR，FKUR] = pKurtosis（流量）;PKUR = FKUR（PKUR == MAX（PKUR））;%计算流量累积和范围。csFlow = cumsum (flow.Data);csFlowRange = max (csFlow)分钟(csFlow);%收集单元格数组中的特征和特征值。%增加流量统计量(平均值、方差等)和公共信号%特性(rms, peak2peak等)到单元阵列。ci = {．..'qmean'，意思（flow.data），．..“qVar”var (flow.Data),．..“qSkewness”偏态(flow.Data),．..“qKurtosis”，峰度（流量数据），．..“qPeak2Peak”peak2peak (flow.Data),．..“qCrest”，峰值均方根（流量数据），．..“qRMS”rms (flow.Data),．..'qmad'，mad（流量数据），．..“qCSRange”csFlowRange,．..“fPeak”fPeak,．..“犁”犁,．..“pMid”pMid,．..“菲格”pHigh,．..“pKurtosis”pKur (1)};结束
另请参阅
相关话题
您还可以从以下列表中选择一个网站：
