MATLAB的24脉波整流器Simulink仿真与故障诊断

本博客来源于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:6subplot(2,3,i)plot(input(30*(i-1)+1,:))axis([-inf inf 3 4.5])grid onylabel('电流幅值/A')xlabel('采样点数/个')if i==1title('正常')elseif i==2title('V1故障')elseif i==3title('V2故障')elseif i==4title('V3故障')elseif i==5title('V4故障')elseif i==6title('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:750iiiiS=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:6subplot(2,3,i)bar(input(30*(i-1)+1,:))%     axis([-inf inf 3 4.5])grid onylabel('电流幅值/A')xlabel('采样点数/个')if i==1title('正常')elseif i==2title('V1故障')elseif i==3title('V2故障')elseif i==4title('V3故障')elseif i==5title('V4故障')elseif i==6title('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 < 2error('ELM:Arguments','Not enough input arguments.');
end  %输入参数必须大于等于2个，否则无法进行建模，因为至少要有输出输入参数才行
if nargin < 3N = size(P_train,2); %如果只有输入输出，默认隐含层神经元为样本数
end
if nargin < 4TF = 'sig';%如果只有输出输入，隐含层神经元数量，默认激活函数为sigmoid函数
end
if nargin < 5TYPE = 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  == 1T1  = 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 TFcase '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  == 1temp_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==0train_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 < 6error('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 TFcase '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 == 1temp_Y = zeros(size(Y));for i = 1:size(Y,2)[max_Y,index] = max(Y(:,i));temp_Y(index,i) = 1;endY = vec2ind(temp_Y); 
end

 

 4 结语

更多内容请点击【专栏】获取。您的点赞是我更新MATLAB1000个案例分析的动力。