本博客来源于CSDN机器鱼，未同意任何人转载。
更多内容，欢迎点击本专栏目录，查看更多内容。
目录
0 引言
1 故障数据采集
2 故障特征提取
3 故障诊断分类
4 结语
本博客内容是在MATLAB2023下完成。
0 引言
对于电力电子电路的故障诊断，作为非专业的电路爱好者，搭建电路非常困难，因此常用simulink仿真来建模，如下图，是一个24脉搏整流器的仿真图。Subsystem是三相整流电路。我用的完整模型可以在【这里】下载。
1 故障数据采集
对于故障诊断数据获取，我们subsystem中的二极管断开连接，就模拟出一种故障类型。
运行下列代码，即可完成一种类型故障数据的采集。一共4个subsystem，每个subsystem有6个二极管，各弄断路一次，一共24种故障，加上正常情况，一共需要运行下面代码25次，得到25个mat文件。记得一定要修改最后保存时候的变量。
clear;clc;close all
format compact
%
%每种故障模式均取30组样本
for i=1:30%每仿真一次就采集一次样本
% 设置负载电阻与电感随机变动，模拟各种负载情况
r=500*(1-0.2*rand);%负载电阻
l=0.1*(1-0.2*rand);%负载电感
set_param('zhengliu24/load','Resistance','r');
set_param('zhengliu24/load','Inductance','l');
sim('zhengliu24');
%保存电流与电压
current=current.signals.values;
voltage=voltage.signals.values;
data{i,1}=current;
data{i,2}=voltage;
end
%% 保存结果 ，每次设置完故障类型，需要将下面这两句的变量名称改一下
%防止覆盖掉
data25=data;
save data25 data25再运行下面的代码，将25个mat文件合并在一起，并给每个样本加上对应的标签。
clear
clc
close all
%% 直接复制，懒得写循环
load data1;
load data2;
load data3;
load data4;
load data4;
load data5;
load data6;
load data7;
load data8;
load data9;
load data10;
load data11;
load data12;
load data13;
load data14;
load data15;
load data16;
load data17;
load data18
load data19
load data20
load data21
load data22
load data23
load data24
load data25
% 要是要电压波形，把1全都改成2就行
data=[data1{:,1} data2{:,1} data3{:,1} ...
data4{:,1} data5{:,1} data6{:,1} ...
data7{:,1} data8{:,1} data9{:,1} ...
data10{:,1} data11{:,1} data12{:,1} ...
data13{:,1} data14{:,1} data15{:,1} ...
data16{:,1} data17{:,1} data18{:,1} ...
data19{:,1} data20{:,1} data21{:,1} ...
data22{:,1} data23{:,1} data24{:,1} data25{:,1} ];
data=data';
label0=1:25;
label1=repmat(label0,30,1);
label=reshape(label1,1,750);
input=data;
save data input label
%% 画图
close all
load data
% 画出正常状态与前5种故障的波形
figure
for i=1:6
subplot(2,3,i)
plot(input(30*(i-1)+1,:))
axis([-inf inf 3 4.5])
grid on
ylabel('电流幅值/A')
xlabel('采样点数/个')
if i==1
title('正常')
elseif i==2
title('V1故障')
elseif i==3
title('V2故障')
elseif i==4
title('V3故障')
elseif i==5
title('V4故障')
elseif i==6
title('V5故障')
end
end
2 故障特征提取
采集完故障数据之后，还需要做故障特征提取，再用提取的特征做故障诊断分类器的训练与测试，代码如下，我们采用快速傅里叶提取频谱作为特征量。
% 特征提取-傅里叶频谱
%
clear
clc
close all
%% 1 加载数据
load data
%% 2 信号加噪 加入均值维0 方差为0.001的高斯白噪声
% 以模拟实际电路中的白噪声与热噪声
m=0;
n=0.001;
data=input+m+n*rand(750,800);
N=800;%采样点数
fs=20000;%采样频率
for iiii=1:750
iiii
S=input(iiii,:);
%% FFT变换并计算模值
Y = fft(S,N); %做FFT变换
Ayy = (abs(Y)); %取模
%% 计算幅频特征
Ayy0=Ayy/(N/2); %换算成实际的幅度
Ayy0(1)=Ayy0(1)/2;
tz(iiii,:)=Ayy0(2:N/2);%同样删除第一分量
end
input=tz;
save data_fft input label
%%
close all
load data_fft
figure
% 画出正常状态与前5种故障的波形
for i=1:6
subplot(2,3,i)
bar(input(30*(i-1)+1,:))
% axis([-inf inf 3 4.5])
grid on
ylabel('电流幅值/A')
xlabel('采样点数/个')
if i==1
title('正常')
elseif i==2
title('V1故障')
elseif i==3
title('V2故障')
elseif i==4
title('V3故障')
elseif i==5
title('V4故障')
elseif i==6
title('V5故障')
end
end
3 故障诊断分类
最后我们采用极限学习机做故障诊断分类，分类前采用PCA做一下降维。
clc;clear;close all; format compact
%% 加载数据
load data_fft
%%
%由于特征维数过高，不利于训练网络，因此我们用PCA进行降维
[PCALoadings,PCAScores,PCAVar] = pca(input);%利用PCA进行降维
percent_explained = 100 * PCAVar / sum(PCAVar);
pareto(percent_explained) %pareto图，只显示前95%的累积分布
xlabel('主成分')
ylabel('贡献率(%)')
title('主成分贡献率')
%
input=PCAScores(:,1:4);%利用pca进行降维至3维，从pareto图看，有4维就够了
%%
output=label';
%随机训练样本，预测样本
[m n]=sort(rand(1,750));
m=500;
P_train=input(n(1:m),:)';
T_train=output(n(1:m),:)';
% P_test=input(n(m+1:end),:)';
% T_test=output(n(m+1:end),:)';
P_test=input(n(end),:)';
T_test=output(n(end),:)';
tic
[IW,B,LW,TF,TYPE,T] = elmtrain(P_train,T_train,100,'sig',1);
toc
train_accuracy=sum(T_train==T)/length(T)
figure
stem(T_train,'r*');hold on
plot(T,'bo');
title('训练集')
%% ELM仿真测试
tic
J = elmpredict(P_test,IW,B,LW,TF,TYPE);
toc
J1=T_test;
disp('展示elm的测试集分类')
test_accuracy=sum(J==J1)/length(J)
figure
stem(J,'r*');hold on
plot(J1,'bo');
title('测试集')function [IW,B,LW,TF,TYPE,T] = elmtrain(P_train,T,N,TF,TYPE)
%P_train=train;T=train_label;N=200;TF='sig';TYPE=1;
% P - Input Matrix of Training Set (R*Q)
% T - Output Matrix of Training Set (S*Q)
% N - Number of Hidden Neurons (default = Q)
% TF - Transfer Function:'sig' for Sigmoidal function (default)',sin' for Sine function,'hardlim' for Hardlim function
% TYPE - Regression (0,default) or Classification (1)
% IW - Input Weight Matrix (N*R)
% B - Bias Matrix (N*1)
% LW - Layer Weight Matrix (N*S)
%回归
% [IW,B,LW,TF,TYPE] = elmtrain(P,T,20,'sig',0)
% Y = elmtrain(P,IW,B,LW,TF,TYPE)
% 分类
% [IW,B,LW,TF,TYPE] = elmtrain(P,T,20,'sig',1)
% Y = elmtrain(P,IW,B,LW,TF,TYPE)
% See also ELMPREDICT
%显然在这个程序中elmtrain跟了6个参数，正常情况下nargin=6，但是有些省略的时候，就要用到下面几句程序进行默认赋值
if nargin < 2
error('ELM:Arguments','Not enough input arguments.');
end %输入参数必须大于等于2个，否则无法进行建模，因为至少要有输出输入参数才行
if nargin < 3
N = size(P_train,2); %如果只有输入输出，默认隐含层神经元为样本数
end
if nargin < 4
TF = 'sig';%如果只有输出输入，隐含层神经元数量，默认激活函数为sigmoid函数
end
if nargin < 5
TYPE = 0;%如果没有定义函数的作用，默认为回归拟合
end
% 其中基本上都是在主程序中定义了那几个参数的，所以上面这几句一般用不到。。。。。。。。。
%%%%%%%%%%%%*****************************
if size(P_train,2) ~= size(T,2)
error('ELM:Arguments','The columns of P and T must be same.');
end
%输出样本数量必须与输出样本数量一致。
[R,Q] = size(P_train);%R=2,Q=1900
if TYPE == 1
T1 = ind2vec(T);
end
%如果定义的是分类，就将训练输出转为向量索引 http://blog.csdn.net/u011314012/article/details/51191006
[S,Q] = size(T1); %S=1,Q=1900
% 随机产生输入权重矩阵，1900*2
IW = rand(N,R) * 2 - 1;
% 随机产生隐层偏置 1900*1
B = rand(N,1);
BiasMatrix = repmat(B,1,Q);
% 计算隐层输出H
tempH = IW * P_train + BiasMatrix;
switch TF
case 'sig'
H = 1 ./ (1 + exp(-tempH));
case 'sin'
H = sin(tempH);
case 'hardlim'
H = hardlim(tempH);
end
% 计算隐层到输出层之间的权重
LW = pinv(H') * T1';
TY=(H'*LW)';
if TYPE == 1
temp_Y=zeros(1,size(TY,2));
for n=1:size(TY,2)
[max_Y,index]=max(TY(:,n));
temp_Y(n)=index;
end
Y_train=temp_Y;
train_accuracy=sum(Y_train==T)/length(T);
end
if TYPE==0
train_accuracy=0;
end
end
%pinv与inv都是用来求矩阵的逆矩阵，但是inv是知道存在逆矩阵的情况下用的
%当不清楚矩阵是否存在逆矩阵的情况时，或者根本不存在逆矩阵，就用PINV来伪逆矩阵function Y = elmpredict(P_test,IW,B,LW,TF,TYPE)
% P - Input Matrix of testing Set (R*Q)
% IW - Input Weight Matrix (N*R)
% B - Bias Matrix (N*1)
% LW - Layer Weight Matrix (N*S)
% TF - Transfer Function:
% 'sig' for Sigmoidal function (default)
% 'sin' for Sine function
% 'hardlim' for Hardlim function
% TYPE - Regression (0,default) or Classification (1)
% Output
% Y - Simulate Output Matrix (S*Q)
% Regression:
% [IW,B,LW,TF,TYPE] = elmtrain(P,T,20,'sig',0)
% Y = elmpredict(P,IW,B,LW,TF,TYPE)
% Classification
% [IW,B,LW,TF,TYPE] = elmtrain(P,T,20,'sig',1)
% Y = elmpredict(P,IW,B,LW,TF,TYPE)
if nargin < 6
error('ELM:Arguments','Not enough input arguments.');
end
% Calculate the Layer Output Matrix H
Q = size(P_test,2);
BiasMatrix = repmat(B,1,Q);
tempH = IW * P_test + BiasMatrix;
switch TF
case 'sig'
H = 1 ./ (1 + exp(-tempH));
case 'sin'
H = sin(tempH);
case 'hardlim'
H = hardlim(tempH);
end
% Calculate the Simulate Output
Y = (H' * LW)';
if TYPE == 1
temp_Y = zeros(size(Y));
for i = 1:size(Y,2)
[max_Y,index] = max(Y(:,i));
temp_Y(index,i) = 1;
end
Y = vec2ind(temp_Y);
end
4 结语
更多内容请点击【专栏】获取。您的点赞是我更新MATLAB1000个案例分析的动力。
