%BP神经网络把这三个传感器的信号特征提取，进行融合，再得出电机状态
clc
clear
close all
%% 第一步 读取数据
feature_index = [1,3,13,14,17];
x=xlsread("data.xlsx");
y=xlsread("10号风机标签.csv");
error_index = find(y~=1);
train_index = [error_index(1),error_index(2),error_index(4),error_index(5),error_index(7),error_index(8)]';
test_index = [error_index(3),error_index(6),error_index(9)]';
y(error_index)
%% 第二步 设置训练数据和预测数据
[m n]=size(x);
input_train =[x(1:20000,feature_index);x(train_index,feature_index)]';
output_train =[y(1:20000,:);y(train_index,:)]';
num = 50;
input_test =[x(20000:20000+num,feature_index);x(error_index,feature_index)]';
output_test =[y(20000:20000+num,:);y(error_index,:)]';
%节点个数
rand('seed',1);%固定随机数，不然每次计算结果有差异
inputnum=size(input_train,1);
hiddennum=[10,4];%隐含层节点
%outputnum=size(output,1);
%% 第三步 训练样本数据归一化
[inputn,inputps]=mapminmax(input_train,0,1);%归一化到[-1,1]之间，inputps用来作下一次同样的归一化
outputn=output_train;
%% 第四步 构建BP神经网络
net=newff(inputn,outputn,hiddennum,{'tansig','purelin'},'trainlm');% 建立模型，传递函数使用purelin，采用梯度下降法训练
W1= net. iw{1, 1};%输入层到中间层的权值
B1 = net.b{1};%中间各层神经元阈值
W2 = net.lw{2,1};%中间层到输出层的权值
B2 = net. b{2};%输出层各神经元阈值
%% 第五步 网络参数配置（ 训练次数，学习速率，训练目标最小误差等）
net.trainParam.epochs=200; % 训练次数
net.trainParam.lr=0.01; % 学习速率，这里设置为0.01
net.trainParam.goal=0.000001; % 训练目标最小误差，这里设置为0.00001
net.trainParam.show = 15;
%% 第六步 BP神经网络训练
[net,tr]=train(net,inputn,outputn);%开始训练，其中inputn,outputn分别为输入输出样本
%% 第七步 测试样本归一化
inputn_test=mapminmax('apply',input_test,inputps);% 对样本数据进行归一化
%% 第八步 BP神经网络预测
test_simu=sim(net,inputn_test);%用训练好的模型进行仿真
%% 第九步 预测结果反归一化与误差计算
%结果整数化
test_simu1=round(test_simu);
test_simu1(test_simu1<1)=1;
test_simu1(test_simu1>4)=4;
error=test_simu1-output_test; %预测值和真实值的误差
acc=length(find(test_simu1==output_test))/length(test_simu1);
%%第十步 真实值与预测值误差比较
figure(1)
stem(output_test,'bo')
hold on
stem(test_simu1,'r*')
hold on
legend('实际故障','预测故障')
xlabel('测试样本')
set(gca,'YTick',0:4)
set(gca,'YTickLabel',{'0' '1-正常' '2-变桨控制故障' '3-叶片故障' '4-CT2故障' })
title("BP-真实值VS预测值");
savefig("1-BP-真实值VS预测值");
[c,l]=size(output_test);
MAE1=sum(abs(error))/l;
MSE1=error*error'/l;
RMSE1=MSE1^(1/2);
disp(['-----------------------误差计算--------------------------'])
disp(['隐含层节点数为',num2str(hiddennum),'时的误差结果如下：'])
disp(['平均绝对误差MAE为：',num2str(MAE1)])
disp(['均方误差MSE为： ',num2str(MSE1)])
disp(['均方根误差RMSE为： ',num2str(RMSE1)])
disp(['分类正确率为/%： ',num2str(acc*100)])
figure(2)
label = {'1-正常' '2-变桨控制故障' '3-叶片故障' '4-CT2故障'};
mat = cfmatrix(output_test,test_simu1);
% mat = [
% sum((test_simu1(1:201)==output_test(1:201)))/201,0,0,0;
% 0,sum((test_simu1(202:204)==output_test(202:204)))/3,0,0;
% sum((test_simu1(205:207)==4))/3,sum((test_simu1(205:207)==3))/3,sum((test_simu1(205:207)==output_test(205:207)))/3,sum((test_simu1(205:207)==1))/3;
% 0,0,sum((test_simu1(205:207)==1))/3,sum((test_simu1(208:210)==output_test(208:210)))/3;
% ];
maxcolor = [191,54,12]; % 最大值颜色
mincolor = [255,255,255]; % 最小值颜色
% 绘制坐标轴
m = length(mat);
imagesc(1:m,1:m,mat)
xticks(1:m)
xlabel('预测类别','fontsize',10.5)
xticklabels(label)
yticks(1:m)
ylabel('实际类别','fontsize',10.5)
yticklabels(label)
% 构造渐变色
mymap = [linspace(mincolor(1)/255,maxcolor(1)/255,64)',...
linspace(mincolor(2)/255,maxcolor(2)/255,64)',...
linspace(mincolor(3)/255,maxcolor(3)/255,64)'];
colormap(mymap)
colorbar()
% 色块填充数字
for i = 1:m
for j = 1:m
text(i,j,num2str(mat(j,i)),...
'horizontalAlignment','center',...
'verticalAlignment','middle',...
'fontname','Times New Roman',...
'fontsize',10);
end
end
% 图像坐标轴等宽
ax = gca;
ax.FontName = 'SimHei';
set(gca,'box','on','xlim',[0.5,m+0.5],'ylim',[0.5,m+0.5]);
set(0,'defaultTextFontName', 'TimesSimSun'); %文字
axis square
savefig("3-BP-混淆矩阵")
