Study on AC load flow using Gauss-seidel method.

Experiment No.: – 10

Study on AC load flow using Gauss-seidel method.

Computer System with MATLAB installed.

Diagram: –

\[\color\magenta{Busdata.m \;file}\] \[\;\] \[\;\] \[\color\green{\% \;Bus\; data\; for\; Load \;Flow \;Analysis.}\] \[\color\red{function\; busdata\; = \;busdata6()}\] \[\color\green{\% \;Returns\; busdata.}\] \[\color\green{\% \; |Bus\; | Type\; | Vsp\; | theta\; | PGi\; | QGi\; | PLi\; | QLi\; | Qmin \;| Qmax\; |}\] \[\color\red{busdata\; =\; [1 \;\;1 \;\;1.05\;\; 0 \;\; 0.0 \;\; 0 \;\; 0 \;\; 0 \;\; 0 \;\; 0;\;\;}\] \[\color\red{2\;\; 2 \;\;1.05\;\;0 \;\;0.5 \;\;0\;\;0\;\;0 \;\;-0.5\;\;1.0;\;\;}\] \[\color\red{3\;\; 2 \;\;1.07 \;\;0 \;\;0.6 \;\;0\;\;0 \;\;0\;\;-0.5\;\;1.5;}\] \[\color\red{4\;\;3\;\;1.0\;\; 0 \;\;0.0\;\;0\;\;0.7 \;\;0.7 \;\; 0 \;\;0;}\] \[\color\red{5\;\;3\;\;1.0\;\;0\;\; 0.0 \;\; 0 \;\; 0.7 \;\;0.7\;\; 0 \;\; 0;}\] \[\color\red{6\;\;3\;\; 1.0 \;\; 0 \;\; 0.0 \;\; 0\;\;0.7\;\;0.7\;\; 0 \;\; 0 ];}\]

% Line Data for Y-Bus Formation.

function linedata = linedata6() % Returns linedata.

% | From | To | R | X | B/2 |
% | Bus | Bus | | | |
linedata = [ 1 2 0.10 0.20 0.02;
1 4 0.05 0.20 0.02;
1 5 0.08 0.30 0.03;
2 3 0.05 0.25 0.03;
2 4 0.05 0.10 0.01;
2 5 0.10 0.30 0.02;
2 6 0.07 0.20 0.025;
3 5 0.12 0.26 0.025;
3 6 0.02 0.10 0.01;
4 5 0.20 0.40 0.04;
5 6 0.10 0.30 0.03;];

y bus:

% Program to form Admittance And Impedance Bus Formation….

function ybus = ybusppg(); % Returns ybus

linedata = linedata6(); % Calling “linedata6.m” for Line Data…
fb = linedata(:,1); % From bus number…
tb = linedata(:,2); % To bus number…
r = linedata(:,3); % Resistance, R…
x = linedata(:,4); % Reactance, X…
b = linedata(:,5); % Ground Admittance, B/2…
z = r + i*x; % Z matrix…

y = 1./z; % To get inverse of each element…

b = i*b; % Make B imaginary…

nbus = max(max(fb),max(tb)); % no. of buses…
nbranch = length(fb); % no. of branches…
ybus = zeros(nbus,nbus); % Initialise YBus…

% Formation of the Off Diagonal Elements…
for k=1:nbranch
ybus(fb(k),tb(k)) = -y(k);
ybus(tb(k),fb(k)) = ybus(fb(k),tb(k));
end

% Formation of Diagonal Elements….
for m=1:nbus
for n=1:nbranch
if fb(n) == m | tb(n) == m
ybus(m,m) = ybus(m,m) + y(n) + b(n);
end
end
end
ybus; % Bus Admittance Matrix
zbus = inv(ybus); % Bus Impedance Matrix
disp(y)

Polar to rectangle:

% Polar to Rectangular Conversion

function rect = pol2rect(r,o) % r = magnitude, o = angle in radians.

rect = rcos(o) + jrsin(o); % rect = real + jimag

Gauss-seidel: –

% Program for Gauss – Seidel Load Flow Analysis

% Assumption, Bus 1 is considered as Slack bus.
clc;
ybus = ybusppg(); % Calling program “ybusppg.m” to get Y-Bus.
busdata = busdata6(); % Calling “busdata6.m” for bus data.

bus = busdata(:,1); % Bus number.
type = busdata(:,2); % Type of Bus 1-Slack, 2-PV, 3-PQ.
V = busdata(:,3); % Initial Bus Voltages.
th = busdata(:,4); % Initial Bus Voltage Angles.
GenMW = busdata(:,5); % PGi, Real Power injected into the buses.
GenMVAR = busdata(:,6); % QGi, Reactive Power injected into the buses.
LoadMW = busdata(:,7); % PLi, Real Power Drawn from the buses.
LoadMVAR = busdata(:,8); % QLi, Reactive Power Drawn from the buses.
Qmin = busdata(:,9); % Minimum Reactive Power Limit
Qmax = busdata(:,10); % Maximum Reactive Power Limit
nbus = max(bus); % To get no. of buses
P = GenMW – LoadMW; % Pi = PGi – PLi, Real Power at the buses.
Q = GenMVAR – LoadMVAR; % Qi = QGi – QLi, Reactive Power at the buses.
Vprev = V;
toler = 1; % Tolerence.
iteration = 1; % iteration starting
while (toler > 0.00001) % Start of while loop
for i = 2:nbus
sumyv = 0;
for k = 1:nbus
if i ~= k
sumyv = sumyv + ybus(i,k)* V(k); % Vk * Yik
end
end
if type(i) == 2 % Computing Qi for PV bus
Q(i) = -imag(conj(V(i))(sumyv + ybus(i,i)V(i)));
if (Q(i) > Qmax(i)) || (Q(i) < Qmin(i)) % Checking for Qi Violation.
if Q(i) < Qmin(i) % Whether violated the lower limit.
Q(i) = Qmin(i);
else % No, violated the upper limit.
Q(i) = Qmax(i);
end
type(i) = 3; % If Violated, change PV bus to PQ bus.
end
end
V(i) = (1/ybus(i,i))((P(i)-jQ(i))/conj(V(i)) – sumyv); % Compute Bus Voltages.
if type(i) == 2 % For PV Buses, Voltage Magnitude remains same, but Angle changes.
V(i) = pol2rect(abs(Vprev(i)), angle(V(i)));
end
end
iteration = iteration + 1; % Increment iteration count.
toler = max(abs(abs(V) – abs(Vprev))); % Calculate tolerance.
Vprev = V; % Vprev is required for next iteration, V(i) = pol2rect(abs(Vprev(i)), angle(V(i)));
end % End of while loop / Iteration

iteration % Total iterations.
V % Bus Voltages in Complex form.
Vmag = abs(V) % Final Bus Voltages.
Ang = 180/pi*angle(V) % Final Bus Voltage Angles in Degree.