LABORATORY MANUAL
POWER SYSTEMS ANALYSIS LABORATORY
EXERCISES
DEPARTMENT OF ELECTRICAL & ELECTRONICS ENGINEERING
Y BUS FORMATION USING MATLAB
AIM: To obtain the Ybus matrix for the given power system using Direct inspection method
and to verify the same using MATLAB.
APPARATUS REQUIRED: Personal Computer with MATLAB software.
FORMATION OF Ybus:
Find the Ybus matrix for the given power system data using Direct inspection method
Sending end Receiving end Reactance values
in ohms
1 2 j0.15
2 3 j0.10
1 3 j0.20
1 4 j0.10
4 3 j0.15
PROGRAM FOR YBUS FORMATION:
% ***** Matlab code for Ybus formation by direct inspection method********** %
clc;
close all;
clear all;
zdata=[1 2 0 0.15;
1 3 0 0.2;
1 4 0 0.1;
2 3 0 0.1;
3 4 0 0.15];
nl=zdata(:,1); % Sending end bus
nr=zdata(:,2); % Receiving end bus
R=zdata(:,3); % Resistance
X=zdata(:,4); % Reactance
nbr=length(zdata(:,1));
nbus = max(max(nl), max(nr));
Z = R + j*X; %branch impedance
y= ones(nbr,1)./Z; %branch admittance
Ybus=zeros(nbus,nbus); % initialize Ybus to zero for k = 1:nbr;
% formation of the off diagonal elements if nl(k) > 0 & nr(k) > 0
Ybus(nl(k),nr(k)) = Ybus(nl(k),nr(k)) - y(k);
Ybus(nr(k),nl(k)) = Ybus(nl(k),nr(k));
end
end
for n = 1:nbus % formation of the diagonal elements for k = 1:nbr
if nl(k) == n | nr(k) == n
Ybus(n,n) = Ybus(n,n) + y(k);
else
end
end
end
Ybus
Z BUS FORMATION USING MATLAB
AIM: To obtain the Zbus matrix for the given power system using Zbus building
algorithm and to verify the same using MATLAB.
APPARATUS REQUIRED: Personal Computer with MATLAB software.
FORMATION OF Zbus:
Find the bus impedance matrix using Zbus building algorithm for the given power system
whose reactance values are as follows.
Sending end Receiving end Reactance values
in ohms
1 0 J0.25
2 1 J0.1
3 1 J0.1
2 0 J0.25
2 3 J0.1
PROGRAM FOR ZBUS FORMATION:
%%%%% final year EEE %%%%%% %%%%%
power system simulation lab%%%% %%%% zbus
building algorithm%%%%%
clc; % clear command window
Clear all; % clear workspace
Close all; % close all windows
%%% Sending bus receving bus reactance value
z=[ 1 0 0.25
2 1 0.1
3 1 0.1
2 0 0.25
2 3 0.1]; % input data
[m n]=size(z); % selecting size of rows and coloums of z data
zbus=[]; % intailizing of zbus with empty matrix
currentbusno=0; % intailizing currentbusno to zero
for a=1:m
[rows cols]=size(zbus); % selecting size of zbus
fb=z(a,1); % selecting sending end bus
tb=z(a,2); % selecting receving end bus
value=z(a,3); % selecting reactance value
newbus=max(fb,tb); % selecting new bus value
ref=min(fb,tb); % selecting ref value
%%%% TYPE-1 MODIFICATION
if newbus>currentbusno && ref==0
zbus=[zbus zeros(rows,1); zeros(1,cols) value];
currentbusno=newbus;
continue
end
%%%% TYPE-2 MODIFICATION
if newbus>currentbusno && ref~=0
zbus=[zbus zbus(:,ref); zbus(ref,:) value+zbus(ref,ref)];
currentbusno=newbus;
continue
end
%%%% TYPE-3 MODIFICATION
if newbus<=currentbusno & ref==0
zbus=zbus-1/(zbus(newbus,newbus)+value)*zbus(:,newbus)*zbus(newbus,:); continue
end
%%%% TYPE-4 MODIFICATION
if newbus<=currentbusno & ref~=0
zbus=zbus-1/(value+zbus(fb,tb)+zbus(tb,tb)-2*zbus(fb,tb))*((zbus(:,fb)-
zbus(:,tb))*(zbus(fb,:)-zbus(tb,:)));
continue
end
end
zbus
GAUSS-SEIDEL LOAD FLOW ANALYSIS USING MATLAB
Aim:- To solve power flow problems by the method of Gauss-Seidel using
MATLAB.
Apparatus Required: -
S.No Apparatus Quantity
1 Personal Computer 1
2 Keyboard 1
3 Mouse 1
4 MATLAB Software 1
Procedure: -
1. Turn on your personal computer.
2. Click on the MATLAB icon of your personal computer.
3. Click the file button and select the new Blank M-file.
4. Type the program on the new M-file for corresponding bus system.
5. After completion of the program, save and run.
6. Note down the line flow and losses.
7. Close the MATLAB tool and turnoff your pc.
% Program for Gauss - Seidel Load Flow Analysis
% Tripavan BTech IV Year, EEE Dept., Mrcet, India, Email
:pavannutulapati@gmail.com
clc;
close all;
clear all;
% From TO Resistance Reactance
linedata = [ 1 2 0.10 0.20;
1 4 0.05 0.20;
1 5 0.08 0.30;
2 3 0.05 0.25;
2 4 0.05 0.10;
2 5 0.10 0.30;
2 6 0.07 0.20;
3 5 0.12 0.26;
3 6 0.02 0.10;
4 5 0.20 0.40;
5 6 0.10 0.30];
nl=linedata(:,1); % Sending end bus
nr=linedata(:,2); % Receiving end bus
R=linedata(:,3); % Resistance
X=linedata(:,4); % Reactance
nbr=length(linedata(:,1));
nbuses = max(max(nl), max(nr));
Z = R + j*X; %branch impedance
y= ones(nbr,1)./Z; %branch admittance
Ybus=zeros(nbuses,nbuses); % initialize Ybus to zero
for k = 1:nbr; % formation of the off diagonal elements
if nl(k) > 0 & nr(k) > 0
Ybus(nl(k),nr(k)) = Ybus(nl(k),nr(k)) - y(k);
Ybus(nr(k),nl(k)) = Ybus(nl(k),nr(k));
end
end
for n = 1:nbuses % formation of the diagonal elements
for k = 1:nbr
if nl(k) == n | nr(k) == n
Ybus(n,n) = Ybus(n,n) + y(k);
else
end
end
end
% |Bus|Type| Vsp |theta |PGi |QGi | PLi | QLi | Qmin | Qmax |
busdata = [ 1 1 1.05 0 0.0 0 0 0 0 0;
2 2 1.05 0 0.5 0 0 0 -0.5 1.0; 3
2 1.07 0 0.6 0 0 0 -0.5 1.5; 4 3
1.0 0 0.0 0 0.7 0.7 0 0; 5 3
1.0 0 0.0 0 0.7 0.7 0 0; 6 3
1.0 0 0.0 0 0.7 0.7 0 0 ]; 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 buse
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)-j*Q(i))/conj(V(i)) - sumyv); % Compute Bus
Voltages.
if type(i) == 2 % For PV Buses, Voltage Magnitude remains same, but Angle
changes.
r=abs(Vprev(i));
o=angle(V(i));
V(i) = r*cos(o) + j*r*sin(o); % rect = real + j*imag
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
disp('***************************************** ')
disp(' Gauss Seidel Load-Flow Study ')
disp('***************************************** ')
disp(' Bus Voltage Angle ')
disp(' no Volts Degree ')
ywz=[ bus(:,1) Vmag Ang ];
disp(ywz)
disp('***************************************** ')
Output:
FAST DECOUPLED LOAD FLOW ANALYSIS USING MATLAB
Aim:- To solve power flow problems by the method of fast decoupled using
MATLAB.
Apparatus Required:-
S.No Apparatus Quantity
1 Personal Computer 1
2 Keyboard 1
3 Mouse 1
4 MATLAB Software 1
Procedure: -
1. Turn on your personal computer.
2. Click on the MATLAB icon of your personal computer.
3. Click the file button and select the new Blank M-file.
4. Type the program on the new M-file for corresponding bus system.
5. After completion of the program, save and run.
6. Note down the line flow and losses.
7. % Program for Fast decoupled Load Flow Analysis
8. clc;
9. close all;
10. clear all;
11.
12. % from to Resistance Reactance%
13. line = [1 4 0.0 0.0576;
14. 4 5 0.017 0.092;
15. 5 6 0.039 0.17;
16. 3 6 0.0 0.0586;
17. 6 7 0.0119 0.1008;
18. 7 8 0.0085 0.072;
19. 8 2 0.0 0.0625;
20. 8 9 0.032 0.161;
21. 9 4 0.01 0.085];
22. nl=line(:,1); % Sending end bus
23. nr=line(:,2); % Receiving end bus
24. R=line(:,3); % Resistance
25. X=line(:,4); % Reactance
26. nbr=length(line(:,1));
27. nbuses = max(max(nl), max(nr));
28. Z = R + j*X; %branch impedance
29. y= ones(nbr,1)./Z; %branch admittance
30. Ybus=zeros(nbuses,nbuses); % initialize Ybus to zero
31. for k = 1:nbr; % formation of the off diagonal elements
32. if nl(k) > 0 & nr(k) > 0
33. Ybus(nl(k),nr(k)) = Ybus(nl(k),nr(k)) - y(k);
34. Ybus(nr(k),nl(k)) = Ybus(nl(k),nr(k));
35. end
36. end
37. for n = 1:nbuses % formation of the diagonal elements
38. for k = 1:nbr
39. if nl(k) == n | nr(k) == n
40. Ybus(n,n) = Ybus(n,n) + y(k);
41. else
42. end
43. end
44. end
45. b = -imag(Ybus);
46. % |Bus| Vsp |theta |PGi |QGi | PLi | QLi | Qmin | Qmax |type
47. bus = [ 1 1.04 0.00 0.00 0.00 0.00 0.00 0.00
0.00
1;
48. 2 1.02533 0.00 1.63 0.00 0.00 0.00 0.00 0.00
2;
49. 3 1.02536 0.00 0.85 0.00 0.00 0.00 0.00 0.00
2;
50. 4 1.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
3;
51. 5 1.00 0.00 0.00 0.00 0.90 0.30 0.00 0.00
3;
52. 6 1.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
3;
53. 7 1.00 0.00 0.00 0.00 1.00 0.35 0.00 0.00
3;
54. 8 1.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
3;
55. 9 1.00 0.00 0.00 0.00 1.25 0.50 0.00 0.00
3];
56.
57. %formation of b' matrtix
58. b1=zeros(nbuses-1,nbuses-1);
59. for i = 1:nbuses-1
60. for j = 1:nbuses-1
61. b1(i,j) =b(i+1,j+1);
62. end
63. end
64. b1;
65. %formation of b" matrtix
66. %assuming all the load buses are at last
67. b2=zeros(nbuses-3,nbuses-3);
68. for i = 1:nbuses-3
69. for j = 1:nbuses-3
70. b2(i,j) =b(i+3,j+3);
71. end
72. end
73. b2;
74. v = bus(:,2);
75. del = bus(:,3);
76. Pg = bus(:,4);
77. Qg = bus(:,5);
78. Pd = bus(:,6);
79. Qd = bus(:,7);
80. Pspec = Pg-Pd;
81. Qspec = Qg-Qd;
82. iter = 1;
83. slack = 1;
84. tolerance = .01;
85. flag=1;
86. while flag==1;
87.
88. m = real(Ybus);
89. n = imag(Ybus);
90. P = zeros(nbuses,1);
91. Q = zeros(nbuses,1);
92. iter= iter+1;
93.
94.
95. for i=1:nbuses %finding bus real and reactive power
96. for j=1:nbuses
97. P(i) = P(i)+ (v(i)*v(j)*(m(i,j)*cos(del(i)-
del(j))+n(i,j)*sin(del(i)-del(j))));
98. Q(i) = Q(i)+ (v(i)*v(j)*(m(i,j)*sin(del(i)-del(j))-
n(i,j)*cos(del(i)-del(j))));
99. end 100.
end 101.
102. P
103. Q
104. %finding del P by v
105. for i=1:(nbuses-1)
106. if(i<slack)
107. delP(i,1)= Pspec(i)-P(i);
108. else
109. delP(i,1)=(Pspec(i+1)-P(i+1));
110.
111. end
112. delPbyv(i,1)=delP(i,1)/v(i,1);
113. end
114.
115. %finding del Q by v
116. c=0;
117. for i=1:nbuses
118. if bus(i,10)==3
119. c=c+1;
120. delQ(c,1)= (Qspec(i)-Q(i));
121. delQbyv(c,1)= delQ(c,1)/v(i,1);
122. end
123. end
124.
125. if max(abs(delP))>tolerance | max(abs(delQ))>tolerance
126. flag=1; % tolerance check
127. else
128. flag=0;
129. end
130.
131. %calc correction vector
132. deldel = inv(b1)*delPbyv;
133. delv = inv(b2)*delQbyv;
134.
135. %updating values
136. for i=1:(nbuses-1)
137. del(i+1,1)= del(i+1,1)+deldel(i,1);
138. end
139.
140. c=0;
141. for i=1:nbuses
142. if bus(i,10)==3
143. c=c+1;
144. v(i,1)=v(i,1)+delv(c,1);
145. end
146. end
147.
148. iter
149. v
150. del
151. end
152. disp('***************************************** ')
153. disp(' Fast decoupled Load-Flow Study ')
154. disp('***************************************** ')
155. disp(' Bus Voltage Angle ')
156. disp(' no Volts Degree ')
157. ywz=[ bus(:,1) v del ];
158. disp(ywz)
159. disp('***************************************** ')
Output:
