Identifikasi kandungan harmonisa dan perhitungan daya pada beban non linear

Selamat sore,

Beban non linear pada sistem tenaga listrik akan menyebabkan arus terdistorsi/cacat. Identifikasi cacat arus dapat dilakukan dengan menggunakan FFT. Berikut adalah program untuk mengidentifikasi dan menghitung daya beban non linear dari data arus dan data tegangan pada beban.

% Perhitungan harmonisa beban linear & non linear dengan sumber sinus
% maupun tidak sinus sampai orde ke 40 dari frekuensi fundamental, menghitung
% nilai Irms, Iharmonisa, Vrms, Vharmonisa, THDi, THDv, P, Q, D, S,
% faktor daya, DF, DPF, dari data arus dan tegangan sinus
%% =====================================================
clc
clear
% informasi data
waktu=0.02; % waktu pencuplikan
data=200;   % jumlah data sampel
Fs=data/waktu;
Ts=1/Fs;

% data sinyal
f=50;   % frekuensi sistem dalam Hz
T=1/f;
w=2*pi*f;
t= (0:data-1)*Ts;
% sinyal tegangan
Vmax0=0;    % tegangan DC
vt0=Vmax0;
Vmax1=100;  % tegangan frekuensi dasar
sudV1=0;
vt1=Vmax1*cos(w*t+sudV1/180*pi);
Vmax2=0;  % tegangan 2 x frekuensi dasar
sudV2=0;
vt2=Vmax2*cos(w*t+sudV2/180*pi);
vt=vt0+vt1+vt2;

% sinyal arus
Imax0=8;    % arus DC
it0=Imax0;
Imax1=15;    % arus frek fundamental
sudI1=-30;    % dalam derajat
it1=Imax1*cos(w*t+sudI1/180*pi);
Imax2=6;    % arus 2 x frek fundamental
sudI2=-45;    % dalam derajat
it2=Imax2*cos(2*w*t+sudI2/180*pi);
Imax3=2;    % arus 3 x frek fundamental
sudI3=-60;    % dalam derajat
it3=Imax3*cos(3*w*t+sudI3/180*pi);
it=it0+it1+it2+it3;

% tampilkan data pada grafik
figure(1);
subplot(2,1,1), plot(t,vt), grid, title('Sinyla Tegangan'),
xlabel('waktu (s)'); ylabel('tegangan (V)');
subplot(2,1,2), plot(t,it), grid, title('Sinyla Arus'),
xlabel('waktu (s)'); ylabel('Arus (Amp.)');

%% Hitung nilai RMS dari sinyal
% hitung rms dengan fungsi inttrazku.m
% sinyal tegangan
Vrmshit=sqrt(1/T*inttrazku(t,vt.^2));
HasilV={'Vrms hit' Vrmshit };
Irmshit=sqrt(1/T*inttrazku(t,it.^2));
HasilI={'Irms hit'; Irmshit};

%% =======================================
% Proses identifikasi / konstruksi sinyal dari data
% FFT untuk mendapatkan magnitidu dan sudut dari sinyal
% FFT sinyal tegangan
L=data;
Y = fft(vt);
f = Fs*(0:(L/2))/L;
% magnitude tegangan
P2 = abs(Y/L);
vP1 = P2(1:L/2+1);
vP1(2:end-1) = 2*vP1(2:end-1);
% sudut tegangan
sP2 = angle(Y/L)*180/pi;
svP1 = sP2(1:L/2+1);
svP1(2:end-1) = svP1(2:end-1);
sVasli= svP1;
% menampilkan hasil
figure(2);
subplot(2,1,1);
bar(f(1:40),vP1(1:40)), grid,
title('Single-Sided Amplitude Spectrum of v(t)')
xlabel('f (Hz)')
ylabel('|V(f)|')

% FFT sinyal arus
Y = fft(it);
f = Fs*(0:(L/2))/L;
% magnitude arus
P2 = abs(Y/L);
iP1 = P2(1:L/2+1);
iP1(2:end-1) = 2*iP1(2:end-1);
% sudut arus
sP2 = angle(Y/L)*180/pi;
siP1 = sP2(1:L/2+1);
siP1(2:end-1) = siP1(2:end-1);
sIasli=siP1;
% menampilkan hasil
subplot(2,1,2);
bar(f(1:40),iP1(1:40)), grid,
title('Single-Sided Amplitude Spectrum of i(t)')
xlabel('f (Hz)')
ylabel('|I(f)|')

%% Hasil konstruksi data apakah sesuai dengan data asli
% cek pembacaan fft sama dengan data untuk magnitude dan sudut
Hasil_tegangan=[vP1(1:40); sVasli(1:40)]
Hasil_arus=[iP1(1:40); sIasli(1:40)]

%% Check nilai RMS berdasarkan konstruksi FFT
Vrmsfft=sqrt(vP1(1).^2+sum((vP1(2:40)/sqrt(2)).^2));
vpersenerrorfft=abs(Vrmsfft-Vrmshit)/Vrmsfft*100;
Irmsfft=sqrt(iP1(1).^2+sum((iP1(2:40)/sqrt(2)).^2));
ipersenerrorfft=abs(Irmsfft-Irmshit)/Irmsfft *100;
HasilV={'Vrms hit' 'Vrms fft' 'kesalahan (%)'; ...
    Vrmshit Vrmsfft vpersenerrorfft}
HasilI={'Irms hit' 'Irms fft' 'kesalahan (%)'; ...
    Irmshit Irmsfft ipersenerrorfft}

%% Perhitungan THD, P, Q, D, S, faktor daya, DF, DPF
rmsV1=vP1(2)/sqrt(2);
THDv=sqrt((Vrmsfft^2-rmsV1^2)/rmsV1^2)*100;
rmsI1=iP1(2)/sqrt(2);
THDi=sqrt((Irmsfft^2-rmsI1^2)/rmsI1^2)*100;
P=vP1(1)*iP1(1)+sum((vP1(2:40).*iP1(2:40))/2.*cos((sVasli(2:40)-...
    sIasli(2:40))*pi/180)); % power
S=Vrmsfft*Irmsfft;       % Apparent power
% cek lag/lead
sudut=(sVasli(2)- sIasli(2))*pi/180;
DPF=cos(sudut);
DF=sqrt(1/(1+(THDi/100).^2));  % distotion faktor
PF=DF*DPF;               % displacement power faktor
if sudut < 0,
    jbeb='lead';
elseif sudut == 0,
    jbeb='unity';
else
    jbeb='lag';
end
D=vP1(2)/sqrt(2)*sqrt(Irmsfft^2-rmsI1^2);
Q=vP1(2)*iP1(2)/2*sin(sudut);
if sudut >= 0,
    Q=Q;
else
    Q=-Q;
end
Hasil={'S' 'P' 'Q' 'D' 'PF' 'DF' 'DPF' 'lead/lag' 'THDv (%)' 'THDi (%)';...
    S P Q D PF DF DPF jbeb THDv THDi}

Semoga bermanfaat.

I M Ari N.