基于电能质量分类的 ML 和 DWT（Matlab实现）

目录

1 Matlab代码

2 结果 

---

本文直接分享代码：

1 Matlab代码

clc
clear
n=250;
%--------------------------------------------------------------------------------------------------------
%% 数学模型的生成
%--------------------------------------------------------------------------------------------------------

for i=1:n
    warning off;
    
%% 正常的50 Hz 正弦波
t=[0 :0.0001:0.4];
k=rand(1,1)*range([1 1.5])+min([1 1.5]);
y=k*(sin((2*pi*50)*t));

figure(1)
subplot(2,1,1);
plot(t,y)
title('Pure 50 Hz Sine wave')
xlabel ('Time (sec)');
ylabel ('Amplitude');
hold on
Pure_Sinal=[t,y]';


%% SAG
%alpha ranges 0.1 to 0.9
t=[0 :0.0001:0.4];
alpha=rand(1,1)*range([0.1 0.8])+min([0.1 0.8]);
t1=0.05;
t2=0.15;
k=rand(1,1)*range([1 1.5])+min([1 1.5]);
y=k*(1-alpha*((heaviside(t-t1)-heaviside(t-t2)))).*sin((2*pi*50)*t);

subplot(2,1,2);
plot(t,y);
title('Sag disturbance');
xlabel ('Time (sec)');
ylabel ('Amplitude');
hold on
Sag=[t,y]';

%% SWELL
%alpha ranges 0.1 to 0.9
t=[0 :0.0001:0.4];
alpha=rand(1,1)*range([0.1 0.8])+min([0.1 0.8]);
t1=0.05;
t2=0.15;
k=rand(1,1)*range([1 1.5])+min([1 1.5]);
y=k*(1+ alpha*((heaviside(t-t1)-heaviside(t-t2)))).*sin((2*pi*50)*t);

figure(2)
subplot(2,1,1);
plot(t,y);
title('Swell disturbance');
xlabel ('Time (sec)');
ylabel ('Amplitude');
hold on
Swell=[t,y]';


%暂停
%alpha ranges 0.9 to 1
t=[0 :0.0001:0.4];
alpha=rand(1,1)*range([0.93 1])+min([0.93 1]);
t1=0.05;
t2=0.15;
k=rand(1,1)*range([1 1.5])+min([1 1.5]);
y=k*(1-alpha*((heaviside(t-t1)-heaviside(t-t2)))).*sin((2*pi*50)*t);

subplot(2,1,2);
plot(t,y);
title('Interruption');
xlabel ('Time (sec)');
ylabel ('Amplitude');
hold on
Interruption=[t,y]';

%% 谐波
%alpha3,alpha5, alpha7 range from 0.05 to 0.15
t=[0 :0.0001:0.4];
alpha3=rand(1,1)*range([0.05 0.15])+min([0.05 0.15]);
alpha5=rand(1,1)*range([0.05 0.15])+min([0.05 0.15]);
alpha7=rand(1,1)*range([0.05 0.15])+min([0.05 0.15]);
alpha1= sqrt(1- alpha3^2-alpha5^2-alpha7^2);
k=rand(1,1)*range([1 1.5])+min([1 1.5]);
y= k*(alpha1* sin(314*t)+ alpha3*sin(3*314*t)+ alpha5*sin(5*314*t)+ alpha7*sin(7*314*t));

figure(3)
subplot(2,1,1);
plot(t,y)
title('Harmonics');
xlabel ('Time (sec)');
ylabel ('Amplitude');
hold on
Harmonics=[t,y]';

%% 暂态脉冲
%fn goes from 300 to 900

fn=500;
amp= rand(1,1)*range([4 7])+min([4 7]);
t1=0.151; 
t2=0.150; 
ty= (t1+t2)/2;
t=[0 :0.0001:0.4];
k=rand(1,1)*range([1 1.5])+min([1 1.5]);
y= k*(sin((2*pi*50)*t)+ amp*(heaviside(t-t2)-heaviside(t-t1)).*exp(-t/ty).*sin(2*pi*fn*t));

subplot(2,1,2);
plot(t,y)
title('Impulsive Transient');
xlabel ('Time (sec)');
ylabel ('Amplitude');
hold on
Impulsive_transient=[t,y]';

%% 振荡瞬态
%fn goes from 300 to 900

fn=rand(1,1)*range([300 500])+min([300 500]);
t=[0 :0.0001:0.4];
amp= 1;
t1=0.255;
t2=0.248; 
ty= (t1+t2)/2;
t=[0 :0.0001:0.4];
k=rand(1,1)*range([1 1.5])+min([1 1.5]);
y= k*(sin((2*pi*50)*t)+ amp*(heaviside(t-t2)-heaviside(t-t1)).*exp(-t/ty).*sin((2*pi*fn)*t));

figure(4)
subplot(2,1,1);
plot(t,y)
title('Oscillatory Transient');
xlabel ('Time (sec)');
ylabel ('Amplitude');
hold on
Oscillatory_transient=[t,y]';

%% SAG+HARMONIC 
t=[0 :0.0001:0.4];
alpha=rand(1,1)*range([0.1 0.8])+min([0.1 0.8]);
alpha3=rand(1,1)*range([0.05 0.15])+min([0.05 0.15]);
alpha5=rand(1,1)*range([0.05 0.15])+min([0.05 0.15]);
alpha7=rand(1,1)*range([0.05 0.15])+min([0.05 0.15]);
alpha1= sqrt(1- alpha3^2-alpha5^2-alpha7^2);
k=rand(1,1)*range([1 1.5])+min([1 1.5]);
y=k*((1-alpha*((heaviside(t-0.05)-heaviside(t-0.15)))).*(alpha1* sin(314*t)+ alpha3*sin(3*314*t)+ alpha5*sin(5*314*t)+ alpha7*sin(7*314*t)));

subplot(2,1,2);
plot(t,y);
title('Sag+Harmonics');
xlabel ('Time (sec)');
ylabel ('Amplitude');
hold on
Sag_harmonic=[t,y]';

%% SWELL+HARMONIC
t=[0 :0.0001:0.4];
alpha=rand(1,1)*range([0.1 0.8])+min([0.1 0.8]);
alpha3=rand(1,1)*range([0.05 0.15])+min([0.05 0.15]);
alpha5=rand(1,1)*range([0.05 0.15])+min([0.05 0.15]);
alpha7=rand(1,1)*range([0.05 0.15])+min([0.05 0.15]);
alpha1= sqrt(1-alpha3^2-alpha5^2-alpha7^2);
k=rand(1,1)*range([1 1.5])+min([1 1.5]);
y=k*((1+alpha*((heaviside(t-0.05)-heaviside(t-0.15)))).*(alpha1* sin(314*t)+ alpha3*sin(3*314*t)+ alpha5*sin(5*314*t)+ alpha7*sin(7*314*t)));

figure(5)
subplot(2,1,1);
plot(t,y)
title('Swell+Harmonics');
xlabel ('Time (sec)');
ylabel ('Amplitude');
hold on
Swell_harmonic=[t,y]';

%% FLICKER
%alpha ranges 0.1 to 2
%beta ranges 5 to 10

t=[0 :0.0001:0.4];
alpha=rand(1,1)*range([0.1 0.15])+min([0.1 0.15]);
beta=rand(1,1)*range([5 7.5])+min([5 7.5]);
k=rand(1,1)*range([1 1.5])+min([1 1.5]);
y=k*((1+alpha*sin(beta*314*t)).*sin(314*t));

subplot(2,1,2);
plot(t,y)
title('Flicker');
xlabel ('Time (sec)');
ylabel ('Amplitude');
hold on
Flicker=[t,y]';

%--------------------------------------------------------------------------------------------------------
%% 特征提取
%--------------------------------------------------------------------------------------------------------


[c1,c2] = dwt(Pure_Sinal,'dmey');% coeficientes discretos
en1 = wentropy(Pure_Sinal,'shannon');% entropia do sinal original
en2 = wentropy(c1,'shannon');% entropia dos coeficiente aproximados (approximation coefficients)
en3 = wentropy(c2,'shannon');% entropia dos coeficientes detalhados (detail coefficients)
% coeff = cwt(vwave); %coeficientes
energy = sqrt(sum(abs(c2).^2,2)); %energia dos coeficientes continuos
x1= mean (energy); %m茅dia da energia
percentage = 100*energy/sum(energy);%energia em porcentagem
[maxpercent,maxScaleIDX] = max(percentage);%maxima porcentagem
r=thd(Pure_Sinal);
T1 (i,:)= table(en1,en2,en3,x1,maxpercent,r,cellstr("pure_sinal"), 'VariableNames', {'Entropy_original', 'Entropy_apcf', 'Entropy_dtcf', 'Mean_energy', 'Max_percent', 'THD', 'Event'});

[c1,c2] = dwt(Sag,'dmey');% coeficientes discretos
en1 = wentropy(Sag,'shannon');% entropia do sinal original
en2 = wentropy(c1,'shannon');% entropia dos coeficiente aproximados (approximation coefficients)
en3 = wentropy(c2,'shannon');% entropia dos coeficientes detalhados (detail coefficients)
% coeff = cwt(vwave); %coeficientes
energy = sqrt(sum(abs(c2).^2,2)); %energia dos coeficientes continuos
x1= mean (energy); %m茅dia da energia
percentage = 100*energy/sum(energy);%energia em porcentagem
[maxpercent,maxScaleIDX] = max(percentage);%maxima porcentagem
r=thd(Sag);
T2(i,:)={en1,en2,en3,x1,maxpercent,r,"sag"};

[c1,c2] = dwt(Swell,'dmey');% coeficientes discretos
en1 = wentropy(Swell,'shannon');% entropia do sinal original
en2 = wentropy(c1,'shannon');% entropia dos coeficiente aproximados (approximation coefficients)
en3 = wentropy(c2,'shannon');% entropia dos coeficientes detalhados (detail coefficients)
% coeff = cwt(vwave); %coeficientes
energy = sqrt(sum(abs(c2).^2,2)); %energia dos coeficientes continuos
x1= mean (energy); %m茅dia da energia
percentage = 100*energy/sum(energy);%energia em porcentagem
[maxpercent,maxScaleIDX] = max(percentage);%maxima porcentagem
r=thd(Swell);
T3(i,:)={en1,en2,en3,x1,maxpercent,r,"swell"};

[c1,c2] = dwt(Interruption,'dmey');% coeficientes discretos
en1 = wentropy(Interruption,'shannon');% entropia do sinal original
en2 = wentropy(c1,'shannon');% entropia dos coeficiente aproximados (approximation coefficients)
en3 = wentropy(c2,'shannon');% entropia dos coeficientes detalhados (detail coefficients)
% coeff = cwt(vwave); %coeficientes
energy = sqrt(sum(abs(c2).^2,2)); %energia dos coeficientes continuos
x1= mean (energy); %dia da energia
percentage = 100*energy/sum(energy);%energia em porcentagem
[maxpercent,maxScaleIDX] = max(percentage);%maxima porcentagem
r=thd(Interruption);
T4(i,:)={en1,en2,en3,x1,maxpercent,r,"interruption"};

[c1,c2] = dwt(Harmonics,'dmey');% coeficientes discretos
en1 = wentropy(Harmonics,'shannon');% entropia do sinal original
en2 = wentropy(c1,'shannon');% entropia dos coeficiente aproximados (approximation coefficients)
en3 = wentropy(c2,'shannon');% entropia dos coeficientes detalhados (detail coefficients)
% coeff = cwt(vwave); %coeficientes
energy = sqrt(sum(abs(c2).^2,2)); %energia dos coeficientes continuos
x1= mean (energy); %m茅dia da energia
percentage = 100*energy/sum(energy);%energia em porcentagem
[maxpercent,maxScaleIDX] = max(percentage);%maxima porcentagem
r=thd(Interruption);
T5(i,:)={en1,en2,en3,x1,maxpercent,r,"harmonics"};

[c1,c2] = dwt(Impulsive_transient,'dmey');% coeficientes discretos
en1 = wentropy(Impulsive_transient,'shannon');% entropia do sinal original
en2 = wentropy(c1,'shannon');% entropia dos coeficiente aproximados (approximation coefficients)
en3 = wentropy(c2,'shannon');% entropia dos coeficientes detalhados (detail coefficients)
% coeff = cwt(vwave); %coeficientes
energy = sqrt(sum(abs(c2).^2,2)); %energia dos coeficientes continuos
x1= mean (energy); %dia da energia
percentage = 100*energy/sum(energy);%energia em porcentagem
[maxpercent,maxScaleIDX] = max(percentage);%maxima porcentagem
r=thd(Impulsive_transient);
T6(i,:)={en1,en2,en3,x1,maxpercent,r,"impulsive_transient"};

[c1,c2] = dwt(Oscillatory_transient,'dmey');% coeficientes discretos
en1 = wentropy(Oscillatory_transient,'shannon');% entropia do sinal original
en2 = wentropy(c1,'shannon');% entropia dos coeficiente aproximados (approximation coefficients)
en3 = wentropy(c2,'shannon');% entropia dos coeficientes detalhados (detail coefficients)
% coeff = cwt(vwave); %coeficientes
energy = sqrt(sum(abs(c2).^2,2)); %energia dos coeficientes continuos
x1= mean (energy); %m茅dia da energia
percentage = 100*energy/sum(energy);%energia em porcentagem
[maxpercent,maxScaleIDX] = max(percentage);%maxima porcentagem
r=thd(Oscillatory_transient);
T7(i,:)={en1,en2,en3,x1,maxpercent,r,"oscillatory_transient"};

[c1,c2] = dwt(Sag_harmonic,'dmey');% coeficientes discretos
en1 = wentropy(Sag_harmonic,'shannon');% entropia do sinal original
en2 = wentropy(c1,'shannon');% entropia dos coeficiente aproximados (approximation coefficients)
en3 = wentropy(c2,'shannon');% entropia dos coeficientes detalhados (detail coefficients)
% coeff = cwt(vwave); %coeficientes
energy = sqrt(sum(abs(c2).^2,2)); %energia dos coeficientes continuos
x1= mean (energy); %m茅dia da energia
percentage = 100*energy/sum(energy);%energia em porcentagem
[maxpercent,maxScaleIDX] = max(percentage);%maxima porcentagem
r=thd(Sag_harmonic);
T8(i,:)={en1,en2,en3,x1,maxpercent,r,"sag_harmonic"};

[c1,c2] = dwt(Swell_harmonic,'dmey');% coeficientes discretos
en1 = wentropy(Swell_harmonic,'shannon');% entropia do sinal original
en2 = wentropy(c1,'shannon');% entropia dos coeficiente aproximados (approximation coefficients)
en3 = wentropy(c2,'shannon');% entropia dos coeficientes detalhados (detail coefficients)
% coeff = cwt(vwave); %coeficientes
energy = sqrt(sum(abs(c2).^2,2)); %energia dos coeficientes continuos
x1= mean (energy); %m茅dia da energia
percentage = 100*energy/sum(energy);%energia em porcentagem
[maxpercent,maxScaleIDX] = max(percentage);%maxima porcentagem
r=thd(Swell_harmonic);
T9(i,:)={en1,en2,en3,x1,maxpercent,r,"swell_harmonic"};

[c1,c2] = dwt(Flicker,'dmey');% coeficientes discretos
en1 = wentropy(Flicker,'shannon');% entropia do sinal original
en2 = wentropy(c1,'shannon');% entropia dos coeficiente aproximados (approximation coefficients)
en3 = wentropy(c2,'shannon');% entropia dos coeficientes detalhados (detail coefficients)
% coeff = cwt(vwave); %coeficientes
energy = sqrt(sum(abs(c2).^2,2)); %energia dos coeficientes continuos
x1= mean (energy); %m茅dia da energia
percentage = 100*energy/sum(energy);%energia em porcentagem
[maxpercent,maxScaleIDX] = max(percentage);%maxima porcentagem
r=thd(Flicker);
T10(i,:)={en1,en2,en3,x1,maxpercent,r,"flicker"};
i
end

%--------------------------------------------------------------------------------------------------------
%CONCATENA脟脙O
%--------------------------------------------------------------------------------------------------------

Matrix = [T1;T2;T3;T4;T5;T6;T7;T8;T9;T10];

%--------------------------------------------------------------------------------------------------------
%% 标准化
%--------------------------------------------------------------------------------------------------------

Matriz_norm=Matrix;

for i=1:n*10
    warning off;
Matriz_norm(i,1)={(Matrix{i,1}-min(Matrix{:,1}))/(max(Matrix{:,1})-min(Matrix{:,1}))};
Matriz_norm(i,2)={(Matrix{i,2}-min(Matrix{:,2}))/(max(Matrix{:,2})-min(Matrix{:,2}))};
Matriz_norm(i,3)={(Matrix{i,3}-min(Matrix{:,3}))/(max(Matrix{:,3})-min(Matrix{:,3}))};
Matriz_norm(i,4)={(Matrix{i,4}-min(Matrix{:,4}))/(max(Matrix{:,4})-min(Matrix{:,4}))};
Matriz_norm(i,5)={(Matrix{i,5}-min(Matrix{:,5}))/(max(Matrix{:,5})-min(Matrix{:,5}))};
Matriz_norm(i,6)={(Matrix{i,6}-min(Matrix{:,6}))/(max(Matrix{:,6})-min(Matrix{:,6}))};
i
end

2 结果