基于粒子群优化深度核极限学习机的故障诊断方法

1-算法原理

2-应用

在博客（点这里）里我们讲述了深度极限学习机，今天我们对其进行改进：原始的深度极限学习机顶层采用极限学习机进行分类，我们将其改成核极限学习机，因为应用了核函数，将低维非线性不可分的数据特征映射到高维线性可分，因此能够划分的更加精准。

1-算法原理

其主要步骤如下：

步骤1：利用极限学习机与自动编码器相结合，形成极限学习机-自动编码器(ELM-AE)，结构如图1，将ELM-AE作为无监督学习的基础单元对输入数据进行贪婪逐层无监督预训练，并保存ELM-AE由最小二乘法获得的输出权值矩阵，

步骤2：然后用这些输出层权重来初始化深度核极限学习机(Deep kernel extreme leaning machine，DKELM)，然后结合标签信息对DKELM的输出层权重进行学习，结构图2所示。

问题引出：由于网络超参数太多，主要有：DKELM各隐含层节点数，ELM-AE的L2 正则化系数，顶层KELM的惩罚系数与核参数，且核函数类型不同，核参数的数量也不一样，因此通过手动选择往往无法得到最优的结果，为此采用粒子群算法对DKELM的超参数进行优化，简称PSO-DKELM。

优化原理：以最小化分类错误率为适应度函数，PSO的目的就是找到一组最优超参数，使得网络的分类错误率最低。

2-应用

为验证上述方法的有效性，本文将PSO-DKELM算法用于凯斯西储的轴承数据故障诊断中。

1）数据准备

采用凯斯西储轴承故障诊断数据，驱动端，48k采样频率，0HP负载，共10类数据，如下所示。从每类数据中各取100个样本，每个样本长度为864个采样点，共得到1000个样本，然后按照7:3划分训练集与测试集。

2）数据处理

对上述采集的样本进行FFT特征提取，提取每个样本的频谱幅值作为最终的特征数据，如下图所示。然后归一化，作为DKELM与PSO-DKELM的输入数据，标签采用onehot编码，即第一类为[1 0 0 0 0 0 0 0 0 0]，第二类为[0 1 0 0 0 0 0 0 0 0]。

3)DKELM分类，参数都是手动设置，随便设的，结果如下

%% DKELM
close all;clear;format compact;format short;clc;rng('default')
%%
load data_feature_fft
[train_x,train_ps]=mapminmax(train_X',-1,1);
test_x=mapminmax('apply',test_X',train_ps);P_train = double(train_x)' ;
P_test  = double(test_x)' ;
T_train = double(train_Y);
T_test  = double(test_Y);
%%
%elm-ae的参数
h=[100 ,50];%各隐含层节点数，里面有几个数就是几个隐含层，即有几个ELM-AE参与无监督预训练，数字代表隐含层节点数
% 比如 h=[100,50],假如输入是200，输出是5类，则DKELM网络结构为200-100-50-5。
% 比如 h=[100,50，20],假如输入是200，输出是5类，则DKELM网络结构为200-100-50-20-5。
TF='sig';%ELM-AE激活函数
lambda=inf;%ELM-AEL2正则化系数
%顶层核极限学习的惩罚系数与核参数
kernel='RBF_kernel';%核函数类型  RBF_kernel lin_kernel poly_kernel wav_kernel
% kernelpara中第一个参数是惩罚系数
if strcmp(kernel,'RBF_kernel')kernelpara=[1 1]; %rbf有一个核参数
elseif strcmp(kernel,'lin_kernel')kernelpara=[1]; %线性没有核参数
elseif strcmp(kernel,'poly_kernel')kernelpara=[1,1,1]; %多项式有2个核参数
elseif strcmp(kernel,'wav_kernel')kernelpara=[1,1,1,1]; %小波核有3个核参数
end
%%
delm=dkelmtrain(P_train,T_train,h,lambda,TF,kernelpara,kernel);
T1=dkelmpredict(delm,kernelpara,kernel,P_train,P_train);[~, J]=max(T1,[],2);
[~, J1]=max(T_train,[],2);
train_accuracy=sum(J==J1)/length(J)
figure
stem(J,'bo');hold on
plot(J1,'r*');
title('训练集分类结果')
legend('实际输出','期望输出')
% 测试集;
T2=dkelmpredict(delm,kernelpara,kernel,P_test,P_train);
[~, J2]=max(T2,[],2);
[~, J3]=max(T_test,[],2);
test_accuracy=sum(J2==J3)/length(J2)
figure
stem(J2,'bo');hold on
plot(J3,'r*');
title('测试集分类结果')
legend('实际输出','期望输出')

由于是随便设置的参数，所以结果不好。也说明做参数优化是有必要的。

4）PSO-DKELM分类，采用PSO优化DKELM（各隐含层节点数，正则化系数，顶层kelm的惩罚系数与核参数），如下，现在只需要设置几个隐含层、ELM-AE激活函数与顶层KELM的核函数类型即可。

%% 优化DKELM超参数,包括各隐含层节点数，正则化系数，顶层kelm的惩罚系数与核参数,
close all;clear;format compact;format short;clc;warning off
%%
load data_feature_fft
[train_x,train_ps]=mapminmax(train_X',-1,1);
test_x=mapminmax('apply',test_X',train_ps);
P_train = double(train_x)' ;
P_test  = double(test_x)' ;
T_train = double(train_Y);
T_test  = double(test_Y);
%% 参数设置
layers=3;%隐含层层数
TF='sig';%ELM-AE激活函数
%顶层核极限学习机参数
kernel='RBF_kernel';%核函数类型  RBF_kernel lin_kernel poly_kernel wav_kernel%%
[x,trace,result]=psofordeepkelm(P_train,P_test,T_train,T_test,layers,TF,kernel); %优化DKELM超参数,包括各隐含层节点数，正则化系数，顶层kelm的惩罚系数与核参数
%将优化结果放到h中
figure
plot(trace)
title('适应度曲线')
xlabel('优化代数')
ylabel('适应度值')hidden=x(1:layers) %各隐含层节点数
lambda=x(layers+1) %L2正则化系数
kernelpara=x(layers+2:end) %顶层kelm惩罚系数与核参数%% 重新训练ML-KELM，看训练集误差
rng('default')
delm=dkelmtrain(P_train,T_train,hidden,lambda,TF,kernelpara,kernel);
T1=dkelmpredict(delm,kernelpara,kernel,P_train,P_train);[~, J]=max(T1,[],2);
[~, J1]=max(T_train,[],2);
train_accuracy=sum(J==J1)/length(J)
figure
stem(J,'bo');hold on
plot(J1,'r*');
title('训练集分类结果')
legend('实际输出','期望输出')
% 测试集;
T2=dkelmpredict(delm,kernelpara,kernel,P_test,P_train);
[~, J2]=max(T2,[],2);
[~, J3]=max(T_test,[],2);
test_accuracy=sum(J2==J3)/length(J2)
figure
stem(J2,'bo');hold on
plot(J3,'r*');
title('测试集分类结果')
legend('实际输出','期望输出')

优化后的适应度曲线如下所示。适应度曲线是一条下降的曲线，表示随着pso优化的，能够找到一组更好的超参数，使得DKELM的分类错误率越来越低。

经过10次优化，得到3个隐含层的节点数分别是74、85、83，ELM-AE的L2正则化系数为542.8677，顶层KELM的惩罚系数为414.3284，核参数为3.9927，采用上述最优超参数进行DKELM建模，得到的测试集分类正确率为99.33%。

4）未来的改进方向

DKELM由于采用了核函数进行，因此不能用于大样本的训练（核矩阵大小为 [训练集样本数 x 训练集样本数]）,以后可以改进这些玩意。