机器学习 深度学习 基于Elman神经网络的刀具剩余使用寿命预测（初学者+matlab代码实现）

1.Elman介绍

        Elman神经网络是 J. L. Elman于1990年首先针对语音处理问题而提出来的，是一种典型的局部回归网络( global feed forward local recurrent)。Elman网络可以看作是一个具有局部记忆单元和局部反馈连接的递归神经网络

        它的主要结构是前馈连接, 包括输入层、 隐含层、 输出层, 其连接权可以进行学习修正;反馈连接由一组“结构 ” 单元构成,用来记忆前一时刻的输出值, 其连接权值是固定的。在这种网络中, 除了普通的隐含层外, 还有一个特别的隐含层,称为关联层 (或联系单元层 ) ;该层从隐含层接收反馈信号, 每一个隐含层节点都有一个与之对应的关联层节点连接。关联层的作用是通过联接记忆将上一个时刻的隐层状态连同当前时刻的网络输入一起作为隐层的输入, 相当于状态反馈。隐层的传递函数仍为某种非线性函数, 一般为 Sigmoid函数, 输出层为线性函数, 关联层也为线性函数。

2.代码实现

2.1数据集来源

刀具数据集来源：试验数据来源于美国纽约预测与健康管理学 会（ＰＨＭ）２０１０年高 速 数 控 机 床 刀 具 健 康 预 测 竞 赛的开放数据。

数据集介绍来源于原文链接：https://blog.csdn.net/qq_40587575/article/details/83351960

实验数据获取的形式是：   试验在上述切削条件下重复进行 ６ 次全寿命周期试验。端面铣削材料为正方形， 每次走刀端     面铣的长度为 １０８ｍｍ 且 每 次 走 刀 时 间 相 等 ， 每次走刀后测量刀具的后刀面磨损量。试验监测数据有ｘ、ｙ 、 ｚ 三向     铣削力信号 ， ｘ 、 ｙ 、 ｚ 三向铣削振动信号以及声发射均方根值。         6次的数据集中  3次实验中有测量铣刀的磨损量，其他3次没有测量，作为比赛的测试集。   系统测量的实验条件和实验方式如下所示：

2.2寿命预测方法流程

1、特征选取：提取六组刀具的切削力信号（一个通道）、振动信号（3个通道）和声发射信号（3个通道），共计7个通道数据，分别提取内个通道的27个特征指标。共计27×7个特征指标，分别上述的方式将特征指标构建成特征矩阵。

2、训练集和测试集选取：随机选取5组刀具数据作为训练集，剩下的一组数据作为测试集。每组数据的长度为315，训练集的长度为1575，测试集的长度为315.

3、数据集标签：刀具的剩余使用寿命标签设置为【315、314、、、、1】

4、构建Elman神经网络，后面就是人工炼丹（调参）

2.3代码实现

2.3.1提取csv文件

%% 提取数据

clc

clear

% 振动信号

file_names = dir('*.csv'); %读取所有的.csv文件

tic; % tic和toc配合使用能够返回程序运行时间

bar = waitbar(0,'计算中...');

for j = 1:length(file_names) %逐一读取文件

data{j} = csvread(file_names(j).name); %依次读取csv格式的数据

str=['计算中...',num2str(100*j/length(file_names)),'%']; % 百分比形式显示处理进程,不需要删掉这行代码就行

waitbar(j/length(file_names),bar,str)

end

toc;

 2.3.2提取特征

%% %% 提取时域特征值

filename = data;

f = struct;

level = 3; % 提前设置存储数据结构

m =size(filename,2);

tic; % tic和toc配合使用能够返回程序运行时间

bar = waitbar(0,'计算中...');

for i = 1:m % 读入通道1的轴承振动信号

for j = 1:7

data = filename{1,i}(:,j); % 需要置于数据文件目录下运

[C,L] = wavedec(data,level,'sym6'); % 小波分解降噪

A = wrcoef('a',C,L,'sym6',level);

% 重构尺度3下的近似小波部分

% % 时域特征

N = length(A);

p1 = mean(A); %均值

y = A-p1;

p2 = sqrt(sum(y.^2)/N); %均方根值,又称有效值(!)

p3 = (sum(sqrt(abs(y)))/N).^2; %方根幅值(!)

p4 = sum(abs(y))/N; %绝对平均值

p5 = sum(y.^3)/N; %歪度

p6 = sum(y.^4)/N; %峭度

p7 = sum((y).^2)/N; %方差

p8 = max(A);%最大值

p9 = min(A);%最小值

p10 = p8-p9;%峰峰值

feature_1.p1(i,j) = p1;

feature_1.p2(i,j) = p2;

feature_1.p3(i,j) = p3;

feature_1.p4(i,j) = p4;

feature_1.p5(i,j) = p5;

feature_1.p6(i,j) = p6;

feature_1.p7(i,j) = p7;

feature_1.p8(i,j) = p8;

feature_1.p9(i,j) = p9;

feature_1.p10(i,j) = p10;

%%以上都是有量纲统计量，以下是无量纲统计量

f1 = p2/p4; %波形指标

f2 = p8/p2; %峰值指标

f3 = p8/p4; %脉冲指标

f4 = p8/p3; %裕度指标

f5 = p6/((p2)^4); %峭度指标

feature_1.f1(i,j) = f1;

feature_1.f2(i,j) = f2;

feature_1.f3(i,j) = f3;

feature_1.f4(i,j) = f4;

feature_1.f5(i,j) = f5;

% 频域特征值提取

fs =20480;

n=0:N-1;

freq=n*fs/N;

f=abs(fft(A,N)*2/N);

freq=freq(1:N/2)';

x=f(1:N/2); % 频率幅值

freq=freq(1:N/2); % 频率值

s1=mean(x); %平均频率

s2=sum(freq.*x)/sum(x);%重心频率

s3=sqrt(sum((freq.^2).*x)/sum(x));%频率均方根

s4=sqrt(sum(((freq-s2).^2).*x)/sum(x));%频率标准差

s5 = f(387,1);

feature_1.s1(i,j) = s1;

feature_1.s2(i,j) = s2;

feature_1.s3(i,j) = s3;

feature_1.s4(i,j) = s4;

feature_1.s5(i,j) = s5;

% 小波bao能量值

wpt=wpdec(A,level,'db4'); %进行3层小波包分解

E = wenergy(wpt); % 提取第三层的小波能量值

feature_1.t1(i,j) = E(:,1);

feature_1.t2(i,j) = E(:,2);

feature_1.t3(i,j) = E(:,3);

feature_1.t4(i,j) = E(:,4);

feature_1.t5(i,j) = E(:,5);

feature_1.t6(i,j) = E(:,6);

feature_1.t7(i,j) = E(:,7);

feature_1.t8(i,j) = E(:,8);

% 显示进度

str=['计算中...',num2str(100*i/m),'%']; % 百分比形式显示处理进程,不需要删掉这行代码就行

waitbar(i/m,bar,str);

end

end

toc;

% 保存数据

clearvars -except feature_1 filename signal T %清除除DATA外的所有变量

2.3.3特征矩阵构建

% save('data.mat');

% 数据集构建（列为特征，行为样本数目

for i =1:7

l(:,i) = feature_1.p1(:,i) ;

l(:,i+7) = feature_1.p2(:,i);

l(:,i+14) = feature_1.p3(:,i) ;

l(:,i+21) = feature_1.p4(:,i) ;

l(:,i+28) = feature_1.p5(:,i) ;

l(:,i+35) = feature_1.p6(:,i) ;

l(:,i+42) = feature_1.p7(:,i) ;

l(:,i+49) = feature_1.p8(:,i) ;

l(:,i+56) = feature_1.p9(:,i) ;

l(:,i+63) = feature_1.p10(:,i) ;

l(:,i+70) = feature_1.f1(:,i) ;

l(:,i+77) = feature_1.f2(:,i) ;

l(:,i+84) = feature_1.f3(:,i) ;

l(:,i+91) = feature_1.f4(:,i) ;

l(:,i+98) = feature_1.f5(:,i) ;

l(:,i+105) = feature_1.s1(:,i) ;

l(:,i+112) = feature_1.s2(:,i) ;

l(:,i+119) = feature_1.s3(:,i) ;

l(:,i+126) = feature_1.s4(:,i) ;

l(:,i+135) = feature_1.t1(:,i) ;

l(:,i+140) = feature_1.t2(:,i) ;

l(:,i+147) = feature_1.t3(:,i) ;

l(:,i+154) = feature_1.t4(:,i) ;

l(:,i+161) = feature_1.t5(:,i) ;

l(:,i+168) = feature_1.t6(:,i) ;

l(:,i+175) = feature_1.t7(:,i) ;

l(:,i+182) = feature_1.t8(:,i) ;

end

%

2.3.4特征筛选

特征筛选思路：将剩余使用寿命标签与特征进行相关性筛选，选取相关性高的特征指标，选取阈值为0.8.

%% 降维

T = (315:-1:1)';

l=l1;

for f= 1:size(l,2)

Ti(:,f) = normalize(l(:,f));

d(1,f) = corr(T,Ti(:,f),'type','Spearman');

end

% 相关性系数

d = (abs(d)); % 将相关系数绝对值

e= sort(d); % 排列特征值的相关性大小

r4=find(d>0.8);

main_feature = d(find(d>=0.8));

% % ±0.80-±1.00 高度相关

M = length(main_feature);

[~,I] = ismember(main_feature,d); % 提取高度相关的特征值

for k = 1:M

p = I(1,k);

value(:,k) = l(:,p); % 提取降维特征

end

2.3.5神经网络实现（读取未经过筛选的特征数据集）

%%数据集降维提取 相关性较高的特征值。

clear

clc

load('D1.mat')

life = (314:-1:0)';

life1 = (314:-1:0)';

life = [life;life1];

life = [life;life1];

life = [life;life1];

life = [life;life1];

for f= 1:size(C1,2)

Ti(:,f) = normalize(C1(:,f));

d(1,f) = corr(life,Ti(:,f),'type','Spearman');

end

% 相关性系数

d = abs(d); % 将相关系数绝对值

[j main_feature] = find(d>=0.8);

% % ±0.80-±1.00 高度相关

M = length(main_feature);

for k = 1:M

p = main_feature(1,k);

c1(:,k) = C1(:,p); % 提取降维特征

c2(:,k) = C2(:,p);

c3(:,k) = C3(:,p);

c4(:,k) = C4(:,p);

c5(:,k) = C5(:,p);

c6(:,k) = C6(:,p);

L1(:,k) = l1(:,p);

L2(:,k) = l2(:,p);

L3(:,k) = l3(:,p);

L4(:,k) = l4(:,p);

L5(:,k) = l5(:,p);

L6(:,k) = l6(:,p);

end

trainx = c1;

trainy = life;

testx = L1;

testy = life1;

%% 创建Elman神经网络

% 包含15个神经元，训练函数为traingdx

net=elmannet(1:2,15,'traingdx');

% 设置显示级别

net.trainParam.show=1;

% 最大迭代次数为2000次

net.trainParam.epochs=2000;

% 误差容限，达到此误差就可以停止训练

net.trainParam.goal=0.00001;

% 最多验证失败次数

net.trainParam.max_fail=5;

% 对网络进行初始化

net=init(net);

%% 网络训练

%训练数据归一化

[trainx1, st1] = mapminmax(trainx');

[trainy1, st2] = mapminmax(trainy');

% 测试数据做与训练数据相同的归一化操作

testx1 = mapminmax('apply',testx',st1);

testy1 = mapminmax('apply',testy',st2);

% 输入训练样本进行训练

[net,per] = train(net,trainx1,trainy1);

%% 测试。输入归一化后的数据，再对实际输出进行反归一化

% 将训练数据输入网络进行测试

train_ty1 = sim(net, trainx1);

train_ty = mapminmax('reverse', train_ty1, st2);

% 将测试数据输入网络进行测试

test_ty1 = sim(net, testx1);

test_ty = mapminmax('reverse', test_ty1, st2);

%% 显示结果

figure(3)

x=1:length(test_ty);

% 显示真实值

plot(x,testy,'b-');

hold on

% 显示神经网络的输出值

plot(x,test_ty,'r--')

legend('真实值','Elman网络输出值')

title('测试数据的测试结果');

mae1=mean(abs(test_ty - testy'));

rmse1 = sqrt(mean((test_ty - testy').^2));

%% 可视化训练及预测结果(过于细致化)

% 显示训练数据的测试结果

% figure(1)

% x=1:length(train_ty);

%

% % 显示真实值

% plot(x,trainy,'b-');

% hold on

% % 显示神经网络的输出值

% plot(x,train_ty,'r--')

%

% legend('真实值','Elman网络输出值')

% title('训练数据的测试结果');

% 显示残差

% figure(2)

% plot(x, train_ty - trainy)

% title('训练数据测试结果的残差')

%

% % 显示均方误差

% mse1 = mse(train_ty - trainy);

% fprintf(' mse_train = \n %f\n', mse1)

%

% % 显示相对误差

% disp(' 训练数据相对误差：')

% fprintf('%f ', (train_ty - trainy)./trainy );

% fprintf('\n')

% figure(3)

% x=1:length(test_ty);

%

% % 显示真实值

% plot(x,testy,'b-');

% hold on

% % 显示神经网络的输出值

% plot(x,test_ty,'r--')

%

% legend('真实值','Elman网络输出值')

% title('测试数据的测试结果');

%

% % % 显示残差

% % figure(4)

% % plot(x, test_ty - testy)

% % title('测试数据测试结果的残差')

%

% % 显示均方误差

%

% mse2 = mse(test_ty - testy);

% rmse1 = sqrt(mean((test_ty - testy').^2));

% % fprintf(' mse_test = \n %f\n', mse2)

% % % 显示相对误差

% % disp(' 测试数据相对误差：')

% % fprintf('%f ', (test_ty - testy)./testy );

% % fprintf('\n')

%

% % % 保存相对误差

% % test_re = test_ty - testy;

% % train_re = train_ty - trainy;

% %

% % test_error = (test_ty - testy)./testy;

% % train_error = (train_ty - trainy)./trainy;

%

数据集表现形式

数据集的预测结果

精彩链接