2022年长三角高校数学建模竞赛

B题齿轮箱故障诊断

原题再现：

齿轮箱是用于增加输出扭矩或改变电机速度的机械装置，被广泛应用于如汽车、输送机、风机等机械设备中。它由两个或多个齿轮组成，其中一个齿轮由电机驱动。电机的轴连接到齿轮箱的一端，并通过齿轮箱的齿轮内部构件，提供由齿轮比确定的输出扭矩和速度。典型的齿轮箱剖面如图 1 所示。在齿轮箱的运行过程中，可以通过加装加速度传感器采集振动信号来判断齿轮箱是否出现异常。本题旨在通过建立相关数学模型对齿轮箱采集到的振动信号进行分析。
在这里插入图片描述
在本题中，我们通过安装在齿轮箱不同部位的四个加速度传感器，采集了 5种状态下齿轮箱的振动信号，具体数据见附件 1。其中表单 gearbox00 为齿轮箱正常工况下采集到的振动信号；表单 gearbox10 为故障状态 1 下采集到的振动信号；表单 gearbox20 为故障状态 2 下采集到的故障信号；表单 gearbox30 为故障状态 3 下采集到的故障信号；表单 gearbox40 为故障状态 4 下采集到的振动信号。信号的采样频率为 6.4kHz。请利用这些数据，建立数学模型解决以下问题：
1、对齿轮箱各个状态下的振动数据进行分析，研究正常和不同故障状态下振动数据的变化规律及差异，并给出刻画这些差异的关键特征。
2、建立齿轮箱的故障检测模型，对其是否处于故障状态进行检测，并对模型的性能进行评价。
3、建立齿轮箱的故障诊断模型，对其处于何种故障状态进行判断，并对模型的性能进行评价。
4、结合所建立的故障检测和诊断模型对附件 2 中另行采集的 12 组测试数据进行检测和诊断分析，将分析结果填写到下表中（注：测试数据中可能存在除以上 4 种故障之外的故障状态，若存在，则将对应的诊断结果标记为：其它故障），并将此表格放到论文的正文中。

整体求解过程概述(摘要)

在工业生产中，齿轮箱作为传递动力的主要部件，其工作状态的好坏是决定机械系统能否高效工作的重要前提。本文在已有数据集的前提下，利用 matlab 软件绘制直方图和求多元统计量来分析各个状态下的振动数据的变化规律及差异，再通过小波降噪、能量分析等方法对数据进行预处理，并用 matlab 软件采用支持向量机方法、XGBoost（机器学习）神经网络预测方法建立齿轮箱故障检测与诊断诊断模型。
针对问题一，旨在深入探究各传感器在各种工况下振动数据变化规律及其差异特征，因此本文将附件一的数据以不同传感器为分类标准分成 4 大组来处理，首先采用数据可视化技术，借助 matlab 工具绘制 20 张不同传感器在不同工况下的时序图和频谱直方图，由图形初步获得数据的主要分布规律及差异。为更深入细节地分析，再通过 matlab 软件求得六种常用多元统计指标即均值、方差、均方差、峭度、脉冲、裕度因子，并对相同传感器的数据进行对比分析，进一步得到数据变化的规律与差异。特别地由峭度、脉冲、裕度的对比，猜测第四个传感器为最大影响因素。
针对问题二，要求建立并评价齿轮箱的故障检测模型，本质是将 5 种工况数据分为故障与非故障两组，即二分类问题，再去建立模型。故本题采用支持向量机方法（svm），运用凸二次规划的最优算法，寻找一个决策超平面，使得二分类效果最优化。再将数据升维到 12 个指标，利用该测试集数据对所得结果进行检验，发现分类情况权重不对等，导致分类层向权重低的空间内入侵，故继续使用 matlab、SPSS 等工具进行模拟加权优化操作，以得到最优的加权比例以及更加科学准确的检测模型。
针对问题三，要求结合前两题特征差异并通过分析建立齿轮箱故障诊断模型，做到对未知样本进行数据故障的检验判别，是对问题二的延伸。故联想采用 XGBoost 分类模型。为排除噪音对于分析的干扰，首先进行数据预处理，采用小波变换、能量分析的方法提取特征数据。接着将所给样本数据以 7：3 的比例划分为训练集和测试集，并使训练集使用 XGBoost 分类模型经过深入学习算法，得到合理的预测模型，最后将测试集代入所得模型，得到训练集与测试集总体正确性在 70%以上，证明该模型的有效性。最后通过模型回归评估得出结论:在 4 个传感器影响中，第 4 个传感器的影响占绝大部分。
针对问题四，需要对 12 组测试数据进行故障诊断，是对问题二以及问题三的模型的具体应用。由于新增一种状况为其他故障，无法先通过问题二模型判别该状况，故本题主要以改进 XGBoost 神经网络预测模型作为基础来增设正常状态下的集合。再建立置信区间估计模型，通过准确值大小先判别是否其他故障，再对剩余组别用修改模型来判别是正常还是四种故障状况之一。最终得到结果为下表：
在这里插入图片描述

模型假设：

1、在解决问题二和问题三时，仅考虑正常工况和已知故障种类，不考虑未知故障种类；
2、不考虑附件二以外其它变量的影响;
3. 假设题目提供数据均为真实测量数值，不存在记录失误等原因造成的错误数据。
4、某个采样部位受到其他部位的影响。

问题分析：

问题一的分析
针对问题一，要求根据题目所给的五个工作表所包含的振动信号数据探究正常工况与故障情况振动信号的变化规律以及其存在的差异，再利用相关方法刻画出存在的差异的关键特征。该工作表数据是在 6.4KHz 采样频率下，通过安装在齿轮箱不同部位的四个加速度传感器，分别采集到的正常工况与四种不同故障情况下齿轮箱的四个不同部位的振动信号数据，故为更加深入直观地探究正常与故障情况的差异，联想到可将所得图像与数据以不同传感器为分类标准均分成四组来分别对比寻找差异。
问题二的分析
针对问题二，要求建立齿轮箱的故障检测模型，并对所建立模型进行评价。即通过将附件 1 中的五组数据分为故障组与非故障组两组，以此来建立故障检测模型。换言之，该题所要建立的模型即将五组情况分为两种类型的模型，是一个二分类模型，再结合问题一的结果，本文初步认定四个传感器的数据有能用于表述故障情况的效度，故联想到应用支持向量机分类模型（svm），支持向量机，它的基本模型是定义在特征空间上的间隔最大的线性分类器，支持向量机的学习策略就是间隔最大化，可形式化为一个求解凸二次规划的问题，支持向量机的学习算法是求解凸二次规划的最优算法，其主要思想是找到一个超平面，使他尽可能多地将两类数据点正确分开，同时也使得两类数据点距离找到的分类平面足够远。从本质上说，支持向量机方法即是一个二分类器（BinaryClassifier）。最后观察利用所建立模型去检验已知样本数据所得到的准确率值，以评价模型的性能。
问题三的分析
针对问题三，要求结合前两题特征差异并通过分析建立齿轮箱故障检验模型，从而对做到对未知样本进行数据故障的检验判别，即建立故障检验预测模型，再对该模型进行评价，是对问题二的延伸。该题需要分析各个故障之间的差异特征，再由此差异特征来预测未知样本数据所属哪一个故障类型，故联想到 XGBoost 算法，它在分类应用中属于多个弱分离器的叠加而成的强分离器，可以根据所给特征差异直接预测所属类别，较适合该题求解。最后，利用验证已知样本数据来计算准确率，作为对模型评价的依据。
问题四的分析
针对问题四，要求根据所建立的故障模型与诊断模型对附件 2 所给的样本数据进行分析判别 12 个测试的工况类型，再填入表格。表格所判别类型多添加了一种其他故障，因此若先用第二题的诊断模型判断其有无故障再采用题三的预测模型得出答案将无法判断其他故障的因素，基于此联想可改进 XGBoost 预测模型在四种故障基础上加入正常工况形成一个新模型。再在预测结果基础上通过数理统计模型先判别为其他故障还是五种工况之一，再对判别后五种工况的情况利用权重比较筛选为哪一组故障或者正常情况。最后在利用所建立的支持向量机检测模型来验证所的结果的准确性。

模型的建立与求解整体论文缩略图

在这里插入图片描述

全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可

部分程序代码：(代码和文档not free)

clc;clear;
figure(1);exam=xlsread(’ 故障 0.xlsx’);
x=exam(:,1)；
a=exam(:,2);b=exam(:,3);
c=exam(:,4);d=exam(:,5);
subplot(2,2,1);plot(x,a,’b’);title(’ 正常 1’);
subplot(2,2,2);plot(x,b,’b’);title(’ 正常 2’);
subplot(2,2,3);plot(x,c,’b’);title(’ 正常 3’);
subplot(2,2,4);plot(x,d,’b’);title(’ 正常 4’);
figure(2);exam=xlsread(’ 故障 1.xlsx’);
x=exam(:,1);
a=exam(:,2);b=exam(:,3);
c=exam(:,4);d=exam(:,5);
subplot(2,2,1);plot(x,a,’b’);title(’ 故障 1’);
subplot(2,2,2);plot(x,b,’b’);title(’ 故障 1’);
subplot(2,2,3);plot(x,c,’b’);title(’ 故障 1’);
subplot(2,2,4);plot(x,d,’b’);title(’ 故障 1’);
figure(3);exam=xlsread(’ 故障 2.xlsx’);
x=exam(:,1);
a=exam(:,2);b=exam(:,3);
c=exam(:,4);d=exam(:,5);
subplot(2,2,1);plot(x,a,’b’);title(’ 故障 2’);
subplot(2,2,2);plot(x,b,’b’);title(’ 故障 2’);
subplot(2,2,3);plot(x,c,’b’);title(’ 故障 2’);
subplot(2,2,4);plot(x,d,’b’);title(’ 故障 2’);
figure(4);exam=xlsread(’ 故障 3.xlsx’);
x=exam(:,1); a=exam(:,2);b=exam(:,3);
c=exam(:,4);d=exam(:,5);
subplot(2,2,1);plot(x,a,’b’);title(’ 故障 3’);
subplot(2,2,2);plot(x,b,’b’);title(’ 故障 3’);
subplot(2,2,3);plot(x,c,’b’);title(’ 故障 3’);
subplot(2,2,4);plot(x,d,’b’);title(’ 故障 3’);
figure(5);exam=xlsread(’ 故障 4.xlsx’);
x=exam(:,1);
a=exam(:,2);b=exam(:,3);
c=exam(:,4);d=exam(:,5);
subplot(2,2,1);plot(x,a,’b’);title(’ 故障 4’);
subplot(2,2,2);plot(x,b,’b’);title(’ 故障 4’);
subplot(2,2,3);plot(x,c,’b’);title(’ 故障 4’);
subplot(2,2,4);plot(x,d,’b’);title(’ 故障 4’);

clc,clear;
data=xlsread(’ 附件 1.xls’,’gearbox00’);
figure(1);y1=data(:,2);hist(y1,10);title(’ 正常’)
figure(2);y2=data(:,3);hist(y2,10);title(’ 正常’)
figure(3);y3=data(:,4);hist(y3,10);title(’ 正常’)
figure(4);y4=data(:,5);hist(y4,100);title(’ 正常’)
data=xlsread(’ 附件 1.xls’,’gearbox10’);
figure(5);q1=data(:,2);hist(q1,10);title(’ 故障 1’)
figure(6);q2=data(:,3);hist(q2,10);title(’ 故障 1’)
figure(7);q3=data(:,4);hist(q3,10);title(’ 故障 1’)
figure(8);q4=data(:,5);hist(q4,100);title(’ 故障 1’)
data=xlsread(’ 附件 1.xls’,’gearbox20’);
figure(9);a1=data(:,2);hist(a1,10);title(’ 故障 2’)
figure(10);a2=data(:,3);hist(a2,10);title(’ 故障 2’)
figure(11);a3=data(:,4);hist(a3,10);title(’ 故障 2’)
figure(12);a4=data(:,5);hist(a4,100);title(’ 故障 2’)
data=xlsread(’ 附件 1.xls’,’gearbox30’);
figure(13);z1=data(:,2);hist(z1,10);title(’ 故障 3’)
figure(14);z2=data(:,3);hist(z2,10);title(’ 故障 3’)
figure(15);z3=data(:,4);hist(z3,10);title(’ 故障 3’)
figure(16);z4=data(:,5);hist(z4,100);title(’ 故障 3’)
data=xlsread(’ 附件 1.xls’,’gearbox40’);
figure(17);w1=data(:,2);hist(w1,10);title(’ 故障 4’)
figure(18);w2=data(:,3);hist(w2,10);title(’ 故障 4’)
figure(19);w3=data(:,4);hist(w3,10);title(’ 故障 4’)
figure(20);w4=data(:,5);hist(w4,100);title(’ 故障 4’)
第一问的均值、方差、均方差、峭度、脉冲、裕度代码：
clc,clear;
data=xlsread(’ 附件 1.xls’,’gearbox00’);
x1=data(:,2);[mean(x1),var(x1),std(x1),kurtosis(x1),(max(x1)-min(x1))/std(x1),
(max(x1)-min(x1))/mean(sqrt(abs(x1)))2
]x2=data(:,3);[mean(x2),var(x2),std(x2),kurtosis(x2),(max(x2)-min(x2))/std(x2),
(max(x2)-min(x2))/mean(sqrt(abs(x2)))2
]
x3=data(:,4);[mean(x3),var(x3),std(x3),kurtosis(x3),(max(x3)-min(x3))/std(x3),
(max(x3)-min(x3))/mean(sqrt(abs(x3)))2
]
x4=data(:,5);[mean(x4),var(x4),std(x4),kurtosis(x4),(max(x4)-min(x4))/std(x4),
(max(x4)-min(x4))/mean(sqrt(abs(x4)))2
]
data=xlsread(’ 附件 1.xls’,’gearbox10’); %故障 1
y1=data(:,2);[mean(y1),var(y1),std(y1),kurtosis(y1),(max(y1)-min(y1))/std(y1),
(max(y1)-min(y1))/mean(sqrt(abs(y1)))2
]
y2=data(:,3);[mean(y2),var(y2),std(y2),kurtosis(y2),(max(y2)-min(y2))/std(y2),
(max(y2)-min(y2))/mean(sqrt(abs(y2)))2
]
y3=data(:,4);[mean(y3),var(y3),std(y3),kurtosis(y3),(max(y3)-min(y3))/std(y3),
(max(y3)-min(y3))/mean(sqrt(abs(y3)))2
]
y4=data(:,5);[mean(y4),var(y4),std(y4),kurtosis(y4),(max(y4)-min(y4))/std(y4),
(max(y4)-min(y4))/mean(sqrt(abs(y4)))2
]
data=xlsread(’ 附件 1.xls’,’gearbox20’);
q1=data(:,2);[mean(q1),var(q1),std(q1),kurtosis(q1),(max(q1)-min(q1))/std(q1),
(max(q1)-min(q1))/mean(sqrt(abs(q1)))2
]
q2=data(:,3);[mean(q2),var(q2),std(q2),kurtosis(q2),(max(q2)-min(q2))/std(q2),
(max(q2)-min(q2))/mean(sqrt(abs(q2)))2
]
q3=data(:,4);[mean(q3),var(q3),std(q3),kurtosis(q3),(max(q3)-min(q3))/std(q3),
(max(q3)-min(q3))/mean(sqrt(abs(q3)))2
]
q4=data(:,5);[mean(q4),var(q4),std(q4),kurtosis(q4),(max(q4)-min(q4))/std(q4),
(max(q4)-min(q4))/mean(sqrt(abs(q4)))2
]
data=xlsread(’ 附件 1.xls’,’gearbox30’);%故障 3
w1=data(:,2);[mean(w1),var(w1),std(w1),kurtosis(w1),(max(w1)-min(w1))/std(w1),
(max(w1)-min(w1))/mean(sqrt(abs(w1)))2
]
w2=data(:,3);[mean(w2),var(w2),std(w2),kurtosis(w2),(max(w2)-min(w2))/std(w2),
(max(w2)-min(w2))/mean(sqrt(abs(w2)))2
]
w3=data(:,4);[mean(w3),var(w3),std(w3),kurtosis(w3),(max(w3)-min(w3))/std(w3),
(max(w3)-min(w3))/mean(sqrt(abs(w3)))2
]
w4=data(:,5);[mean(w4),var(w4),std(w4),kurtosis(w4),(max(w4)-min(w4))/std(w4),
(max(w4)-min(w4))/mean(sqrt(abs(w4)))2
]
data=xlsread(’ 附件 1.xls’,’gearbox40’);%故障 4
z1=data(:,2);[mean(z1),var(z1),std(z1),kurtosis(z1),(max(z1)-min(z1))/std(z1),
(max(z1)-min(z1))/mean(sqrt(abs(z1)))2
]
z2=data(:,3);[mean(z2),var(z2),std(z2),kurtosis(z2),(max(z2)-min(z2))/std(z2),
(max(z2)-min(z2))/mean(sqrt(abs(z2)))2
]
z3=data(:,4);[mean(z3),var(z3),std(z3),kurtosis(z3),(max(z3)-min(z3))/std(z3),
(max(z3)-min(z3))/mean(sqrt(abs(z3)))2
]
z4=data(:,5);[mean(z4),var(z4),std(z4),kurtosis(z4),(max(z4)-min(z4))/std(z4),
(max(z4)-min(z4))/mean(sqrt(abs(z4)))2
]