clear
clc
disp('                                                                                ')
disp('              ____________________________________________________              ')
disp('                                                                                ')
disp('                      STUDI ALIRAN DAYA PADA SISTEM 3 BUS                       ')
disp('                    DENGAN PENYELESAIAN METODE FAST DECOUPLED                   ')
disp('                            OLEH : ZALDI HARDIYANTO                             ')
disp('                         WWW.ELEKTRONIKABERSAMA.WEB.ID                          ')
disp('              ____________________________________________________              ')
disp('                                                                                ')
% Dasar 100 MVA
epsilon = 0.00001;
x = 1;
% Data impedansi pada saluran :
z12 = 0.02 + j*0.04;
z13 = 0.01 + j*0.03;
z23 = 0.0125 + j*0.025;
% Admitansi pada saluran :
y12 = 1/z12;
y13 = 1/z13;
y23 = 1/z23;
% Admitansi Ybus :
Y11 = (y12 + y13);
Y12 = -y12;
Y13 = -y13;
Y21 = Y12;
Y22 = y12 + y23;
Y23 = -y23;
Y31 = -y13;
Y32 = -y23;
Y33 = y13 + y23;
% Pembentukan sudut :
theta11 = atan2(imag(Y11),real(Y11));
theta12 = atan2(imag(Y12),real(Y12));
theta13 = atan2(imag(Y13),real(Y13));
theta21 = atan2(imag(Y21),real(Y21));
theta22 = atan2(imag(Y22),real(Y22));
theta23 = atan2(imag(Y23),real(Y23));
theta31 = atan2(imag(Y31),real(Y31));
theta32 = atan2(imag(Y32),real(Y32));
theta33 = atan2(imag(Y33),real(Y33));
% Elemen-elemen matriks suseptansi untuk mengevaluasi sudut fase delta_thau2 dan thau3 :
B1(1,1)=imag(Y22);
B1(1,2)=imag(Y23);
B1(2,1)=imag(Y32);
B1(2,2)=imag(Y33);
% Matriks invers dari B' :
b1=B1^-1;
% Harga B" :
B2=imag(Y22);
% Beban dalam per unit :
S2 = -(256.6 + j*110.2)/100;
S3 = -(138.6 + j*45.2)/100;
% Bus 1 sebagai slack bus dengan :
V1 = 1.05 + j*0.0;
% Estimasi tegangan awal untuk :
V2 = 1.0 + j*0.0;
V3 = 1.0 + j*0.0;
% Harga thau1 :
thau1 = 0;
% Estimasi thau awal untuk :
thau2 = 0;
thau3 = 0;
iter = 0;
format short g 
while x >= epsilon
    iter = iter + 1;
    % Harga cos dan sin dari sudut :
    sda = cos(theta11);
    sdb = cos(theta12 - thau1 + thau2);
    sdc = cos(theta13 - thau1 + thau3);
    sd1 = cos(theta21 - thau2 + thau1);
    sd2 = cos(theta22);
    sd3 = cos(theta23 - thau2 + thau3);
    sd4 = cos(theta31 - thau3 + thau1);
    sd5 = cos(theta32 - thau3 + thau2);
    sd6 = cos(theta33);
    sd7 = sin(theta21 - thau2 + thau1);
    sd8 = sin(theta22);
    sd9 = sin(theta23 - thau2 + thau3);
    sd10 = sin(theta31 - thau3 + thau1);
    sd11 = sin(theta32 - thau3 + thau2);
    sd12 = sin(theta33);
    sdd = sin(theta11);
    sde = sin(theta12 - thau1 + thau2);
    sdf = sin(theta13 - thau1 + thau3);
    % Harga daya nyata dan daya reaktif pada bus 2 dan 3 :
    P2 = (abs(V2)*abs(V1)*abs(Y21)*sd1) + (abs(V2^2)*abs(Y22)*sd2) + (abs(V2)*abs(V3)*abs(Y23)*sd3);
    P3 = (abs(V3)*abs(V1)*abs(Y31)*sd4) + (abs(V3)*abs(V2)*abs(Y32)*sd5) + (abs(V3^2)*abs(Y33)*sd6);
    Q2 = -(abs(V2)*abs(V1)*abs(Y21)*sd7) - (abs(V2^2)*abs(Y22)*sd8) - (abs(V2)*abs(V3)*abs(Y23)*sd9);
    Q3 = -(abs(V3)*abs(V1)*abs(Y31)*sd10) - (abs(V3)*abs(V2)*abs(Y32)*sd11) - (abs(V3^2)*abs(Y33)*sd12);
    % Sisa daya (Power residuals) :
    deltaP2 = real(S2) - P2;
    deltaP3 = real(S3) - P3;
    deltaQ2 = imag(S2) - Q2;
    deltaQ3 = imag(S3) - Q3;
    % Perhitungan sudut tegangan pada bus 2 dan 3 :
    thau = -b1*[deltaP2/abs(V2); deltaP3/abs(V3)];
    thau2 = thau2 + thau(1);
    thau3 = thau3 + thau(2);
    % Perhitungan harga tegangan pada bus 2 dan 3:
    V = -(1/B2)*[deltaQ2/abs(V2); deltaQ3/abs(V3)];
    V2 = V2 + V(1);
    E3 = V3 + V(2);
    % Tegangan dalam bentuk bilangan kompleks pada bus 2 dan 3 :
    Vkom2 = abs(V2)*(cos(thau2) + j*sin(thau2));
    Vkom3 = abs(E3)*(cos(thau3) + j*sin(thau3));
    % Arus pada saluran :
    I12 = y12*(V1 - Vkom2);
    I21 = -I12;
    I13 = y13*(V1 - Vkom3);
    I31 = -I13;
    I23 = y23*(Vkom2 - Vkom3);
    I32 = -I23;
    % Konjugat dari arus pada saluran :
    Ik12 = conj(I12);
    Ik21 = conj(I21);
    Ik13 = conj(I13);
    Ik31 = conj(I31);
    Ik23 = conj(I23);
    Ik32 = conj(I32);    
    % Daya nyata dan daya reaktif pada bus 1 :
    P1 = (abs(V1^2)*abs(Y11)*sda) + (abs(V1)*abs(Vkom2)*abs(Y12)*sdb) + (abs(V1) *abs(Vkom3)*abs(Y13)*sdc);
    Q1 = -(abs(V1^2)*abs(Y11)*sdd) - (abs(V1)*abs(Vkom2)*abs(Y12)*sde) - (abs(V1)*abs(Vkom3)*abs(Y13)*sdf);
    % Daya kompleks pada bus 1 :
    S1 = P1+j*Q1;
    % aliran daya dalam bentuk bilangan kompleks pada saluran :
    S12 = V1*Ik12*100;
    S21 = Vkom2*Ik21*100;
    S13 = V1*Ik13*100;
    S31 = Vkom3*Ik31*100;
    S23 = Vkom2*Ik23*100;
    S32 = Vkom3*Ik32*100;
    % Rugi-rugi daya dalam bentuk bilangan kompleks pada saluran :
    SL12 = S12+S21;
    SL13 = S13+S31;
    SL23 = S23+S32;
    % Tegangan dalam sudut derajat:
    Vt1 = 180/pi*thau1;
    Vt2 = 180/pi*thau2;
    Vt3 = 180/pi*thau3;
    fprintf('Iterasi %i', iter)
    disp(' ')
    disp('Tegangan dan sudut fase dalam p.u.')
    disp(['V1 = ', num2str(V1),' ', 'Sudut = ', num2str(Vt1)])
    disp(['V2 = ', num2str(V2),' ', 'Sudut = ', num2str(Vt2)])
    disp(['V3 = ', num2str(V3),' ', 'Sudut = ', num2str(Vt3)])
    disp('Daya Kompleks Bus 1 dalam p.u.')
    disp(['S1 = ', num2str(S1)])
    disp('Aliran Daya dan Rugi-rugi Daya dalam satuan MW dan Mvar')
    disp(['S12 = ', num2str(S12),' ', 'S21 = ', num2str(S21)])
    disp(['S13 = ', num2str(S13),' ', 'S31 = ', num2str(S31)])
    disp(['S23 = ', num2str(S23),' ', 'S32 = ', num2str(S32)])
    disp(['SL12 = ', num2str(SL12),' ', 'SL13 = ', num2str(SL13),' ', 'SL23 = ', num2str(SL23)])
    disp(' ')
    x = abs(E3 - V3);
    V3 = E3;
end