Program Mathlab untuk memformulasikan YBus

buatlah program berbasis MATLAB untuk memformulasikan YBus berbasis format input di bawah ini dan tunjukkan hasilnya!

%   Bus Data Format

%       1   bus number (1 to 29997)

%       2   bus type

%               PQ bus            = 1

%               PV bus            = 2

%               reference bus   = 3

%               isolated bus     = 4

%       3   Pd, real power demand (MW)

%       4   Qd, reactive power demand (MVAR)

%       5   Gs, shunt conductance (MW (demanded?) at V = 1.0 p.u.)

%       6   Bs, shunt susceptance (MVAR (injected?) at V = 1.0 p.u.)

%       7   area number, 1-100

%       8   Vm, voltage magnitude (p.u.)

%       9   Va, voltage angle (degrees)

%       10  baseKV, base voltage (kV)

%       11  zone, loss zone (1-999)

%       12  maxVm, maximum voltage magnitude (p.u.)

%       13  minVm, minimum voltage magnitude (p.u.)

%

%   Generator Data Format

%       1   bus number

%       2   Pg, real power output (MW)

%       3   Qg, reactive power output (MVAR)

%       4   Qmax, maximum reactive power output (MVAR)

%       5   Qmin, minimum reactive power output (MVAR)

%       6   Vg, voltage magnitude setpoint (p.u.)

%       7   mBase, total MVA base of this machine, defaults to baseMVA

%       8   status, 1 – machine in service, 0 – machine out of service

%       9   Pmax, maximum real power output (MW)

%       10  Pmin, minimum real power output (MW)

%

%   Branch Data Format

%       1   f, from bus number

%       2   t, to bus number

%       3   r, resistance (p.u.)

%       4   x, reactance (p.u.)

%       5   b, total line charging susceptance (p.u.)

%       6   rateA, MVA rating A (long term rating)

%       7   rateB, MVA rating B (short term rating)

%       8   rateC, MVA rating C (emergency rating)

%       9   ratio, transformer off nominal turns ratio ( = 0 for lines )

%           (taps at ‘from’ bus, impedance at ‘to’ bus, i.e. ratio = Vf/Vt)

%       10  angle, transformer phase shift angle (degrees)

%       11  initial branch status, 1 – in service, 0 – out of service

%% system MVA base

baseMVA = 100.0000;

%% bus data

%bus type   Pd    Qd     Gs      Bs        area                 Vm       Va          baseKV    zone          Vmax     Vmin

bus = [

1   1    40.0   15.0    0.0     0.0           1              1.0000     0.0000     230.0000   1              1.0500  0.9500;

2   3    25.0   10.0    0.0     0.0           1              1.0000     0.0000     230.0000   1              1.0500  0.9500;

3   2    60.0   30.0    0.0     0.0           1              1.0000     0.0000     230.0000    1              1.0500  0.9500;

4   1    70.0   25.0    0.0     0.0           1              1.0000     0.0000     230.0000   1              1.0500  0.9500;

];

%% generator data

% bus  Pg        Qg      Qmax      Qmin                      Vsp            base     status           Pmax                   Pmin

gen = [

2   85.00  19.22  150.0000 -150.0000                1.050       100.0000    1             200.0000              50.0000;

3  110.00  67.43 150.0000 -150.0000                1.050       100.0000    1             200.0000              50.0000;

];

%% branch data

%         fbus          tbus               r              x                b/2              ratea       rateb      ratec           ratio     angle    status

branch = [

1              2              0.20         0.60         0.0010     40.0000  40.0000   250.0000   0.0000   0.0000     1;

1              2              0.20         0.60         0.0010     40.0000  40.0000   250.0000   0.0000   0.0000     1;

1              3              0.10         0.20         0.0000     80.0000  80.0000   150.0000   0.0000   0.0000     1;

2              3              0.07         0.20         0.0000     60.0000  60.0000   150.0000   0.0000   0.0000     1;

2              4              0.05         0.25         0.0000     100.000   100.000   150.0000   0.0000   0.0000     1;

3              4              0.04         0.20         0.0010     50.0000  50.0000   150.0000   0.0000   0.0000     1;

3              4              0.04         0.20         0.0010     50.0000  50.0000   150.0000   0.0000   0.0000     1;

];

Jalankan program di bawah ini di Mathlab

function[Ybus] = ybus(branch)

%% bus data

%bus type   Pd    Qd     Gs      Bs        area                 Vm       Va          baseKV    zone          Vmax     Vmin

bus = [

1   1   40.0   15.0    0.0  0.0  1  1.0000  0.0000  230.0000    1     1.0500  0.9500;

2   3   25.0   10.0    0.0  0.0  1  1.0000  0.0000  230.0000    1     1.0500  0.9500;

3   2   60.0   30.0    0.0  0.0  1  1.0000  0.0000  230.0000    1     1.0500  0.9500;

4   1   70.0   25.0    0.0  0.0  1  1.0000  0.0000  230.0000    1     1.0500  0.9500;

];

%% generator data

% bus  Pg        Qg      Qmax      Qmin                      Vsp            base     status           Pmax                   Pmin

gen = [

2   85.00  19.22  150.0000 -150.0000  1.050  100.0000     1     200.0000  50.0000;

3  110.00  67.43  150.0000 -150.0000  1.050  100.0000     1     200.0000  50.0000;

];

%% branch data

%fbus tbus   r   x     b/2              ratea       rateb      ratec           ratio     angle    status

branch = [

1  2  0.20  0.60  0.0010   40.0000  40.0000  250.0000   0.0000  0.0000  1;

1  2  0.20  0.60  0.0010   40.0000  40.0000  250.0000   0.0000  0.0000  1;

1  3  0.10  0.20  0.0000   80.0000  80.0000  150.0000   0.0000  0.0000  1;

2  3  0.07  0.20  0.0000   60.0000  60.0000  150.0000   0.0000  0.0000  1;

2  4  0.05  0.25  0.0000   100.000  100.000  150.0000   0.0000  0.0000  1;

3  4  0.04  0.20  0.0010   50.0000  50.0000  150.0000   0.0000  0.0000  1;

3  4  0.04  0.20  0.0010   50.0000  50.0000  150.0000   0.0000  0.0000  1;

];

busfrom=branch(:,1); busto=branch(:,2); Resistansi=branch(:,3); reaktansi=branch(:,4);

nbr=length(branch(:,1)); nbus = max(max(busfrom), max(busto));

%Rumus Z

Z = Resistansi + j*reaktansi;

y= ones(nbr,1)./Z;

Ybus=zeros(nbus,nbus);

for n = 1:nbus  %fungsi perulangan untuk elemen diagonal

for a = 1:nbr

if busfrom(a) == n || busto(a) == n

Ybus(n,n) = Ybus(n,n) + y(a);

else, end

end

end

for a = 1:nbr;   %fungsi perulangan untuk elemen selain diagonal

if busfrom(a) > 0 && busto(a) > 0

Ybus(busfrom(a),busto(a)) = Ybus(busfrom(a),busto(a)) – y(a);

Ybus(busto(a),busfrom(a)) = Ybus(busfrom(a),busto(a));

end

end

Hasilnya setelah dirun

ans =

Columns 1 through 3

3.0000 – 7.0000i        -1.0000 + 3.0000i      -2.0000 + 4.0000i        0

-1.0000 + 3.0000i        3.3283 -11.3005i      -1.5590 + 4.4543i   -0.7692 + 3.8462i

-2.0000 + 4.0000i       -1.5590 + 4.4543i       5.4821 -18.0697i  -1.9231 + 9.6154i

0                                 -0.7692 + 3.8462i      -1.9231 + 9.6154i     2.6923 -13.4615i

Hasilnya sesuai dengan perhitungan manual

Image

Image