clc
clear all
close all
format shortg
%% Bus Data
% "Type": 1-slack bus; 2-PV bus; 3-PQ bus;
% Base 100MVA, 230KV, convert P and Q into per unit
% |Bus | Type | Vsp | theta | PGi | QGi | PLi | QLi |
busdata = [ 1 1 1.00 0 0.0 0 50/100 30.99/100;
2 3 1.00 0 0.0 0 170/100 105.35/100;
3 3 1.00 0 0.0 0 200/100 123.94/100;
4 2 1.02 0 318/100 0 80/100 49.58/100];
%% Line Data
% | From | To | R | X | Y/2 |
% | Bus | Bus | | | |
linedata = [ 1 2 0.01008 0.05040 0.05125;
1 3 0.00744 0.03720 0.03875;
2 4 0.00744 0.03720 0.03875;
3 4 0.01272 0.06360 0.06375];
%% Formulate Admittance and Impedance Matrix
fb = linedata(:,1); % From bus number
tb = linedata(:,2); % To bus number
r = linedata(:,3); % Resistance, R
x = linedata(:,4); % Reactance, X
b = j*linedata(:,5); % Ground Admittance, Y/2...
z = r + j*x; % Z matrix
y = 1./z; % Y=inv(Z)
nbus = max(max(fb),max(tb)); % Number of buses
nbranch = length(fb); % Number of branches...
ybus = zeros(nbus,nbus); % Initialize Y-Bus...
% Number of PV buses
pvbus = 0;
for n=1:nbus
if busdata(n,2)==2
pvbus=pvbus+1;
end
end
% Number of PQ buses
pqbus = 0;
for n=1:nbus
if busdata(n,2)==3
pqbus=pqbus+1;
end
end
% Yjj = Sum of admittances connected to node j
% Yjk = - (Sum of the admittances connected between node j and node k)
% 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
% 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
% Check the bus admittance matrix "ybus" with your solutions
ybus % Bus Admittance Matrix
%% Gauss-Seidel Power Flow
% Bus 1 is a slack bus (reference)
iteration = 1;
type = busdata(:,2);
V = busdata(:,3); % Initial bus voltages
th = busdata(:,4); % Initial bus voltage angles.
GenMW = busdata(:,5); % PGi, Real Power injected into the bus
GenMVAR = busdata(:,6); % QGi, Reactive Power injected into the bus
LoadMW = busdata(:,7); % PLi, Real Power drawn from the bus
LoadMVAR = busdata(:,8); % QLi, Reactive Power drawn from the bus
P = GenMW - LoadMW; % Pi = PGi - PLi
Q = GenMVAR - LoadMVAR; % Qi = QGi - QLi
toler = 1; % Tolerence
IterMax = 100; % Maximum # of iterations
Vpe = V;
Va = V(2:end);
Qa = Q(2:end);
tha = th(2:end);
Pa = P(2:end);
i=0;
for k=1:nbus
if busdata(k,2)== 2
Va(k-1-i) = [];
Qa(k-1-i) = [];
%%%%%%%%%%%%%%%%%%%%%
i=i+1;
%%%%%%%%%%%%%%%%%%%%%
end
end
theta_V = [tha;Va];
for iteration = 1:1:IterMax
% Initialize Jacobian Matrix
J = zeros((2*nbus-2-pvbus),(2*nbus-2-pvbus));
J11=zeros(nbus-1,nbus-1);
J12=zeros(nbus-1,pqbus);
J21=zeros(pqbus,nbus-1);
J22=zeros(pqbus,pqbus);
% Delta P
for k = 2:nbus
YVV = 0;
for n = 1:nbus
YVV = YVV + (abs(ybus(k,n))* abs(V(n))*abs(V(k)))*cos((th(k))-th(n)-angle(ybus(k,n)));
end
Pc(k-1,1) = YVV;
deltaP(k-1,1)=Pa(k-1)-Pc(k-1);
end
Pa
Pc
deltaP
% Delta Q
for k = 2:(nbus-pvbus)
YVq = 0;
if busdata(k,2) ~= 2
for n = 1:nbus
YVq = YVq + (abs(ybus(k,n))*abs(V(n))*abs(V(k)))*sin((th(k))-th(n)-angle(ybus(k,n)));
end
Qc(k-1,1) = YVq;
deltaQ(k-1,1)=Qa(k-1)-Qc(k-1);
end
end
% Delta PQ Matrix
deltaPQ = [deltaP;deltaQ];
% Calculate J11
for k=2:nbus
for n=2:nbus
if k~=n
J11(k-1,n-1)= abs(V(k))*abs(ybus(k,n))*abs(V(n))*sin(-angle(ybus(k,n))+th(k)-th(n));
end
if k==n
for nn=1:nbus
if nn~=k
J11(k-1,k-1) = J11(k-1,k-1)-abs(V(k))*abs(ybus(k,nn))*abs(V(nn))*sin(angle(-ybus(k,nn))+th(k)-th(nn));
end
end
end
end
end
% Calculate J12
for k=2:nbus
for n=2:nbus-pvbus
if busdata(n,2)==3
if k~=n
J12(k-1,n-1)= abs(ybus(k,n))*abs(V(k))*cos(angle(ybus(k,n))+th(n)-th(k));
end
yv12=0;
if k==n
for nn=1:nbus
if nn~=k
yv12 = yv12+abs(ybus(k,nn))*abs(V(nn))*cos(angle(ybus(k,nn)+th(nn)-th(k)));
end
end
J12(k-1,k-1) = 2*abs(V(k))*abs(ybus(k,k))*cos(angle(ybus(k,k)))+yv12;
end
end
end
end
% Calculate J21
for k=2:nbus-pvbus
for n=2:nbus
if busdata(k,2)==3
if k~=n
J21(k-1,n-1)= -abs(V(k))*abs(ybus(k,n))*abs(V(n))*cos(-angle(ybus(k,n))-th(n)+th(k));
end
if k==n
for nn=1:nbus
if nn~=k
J21(k-1,k-1) = J21(k-1,k-1)+abs(V(k))*abs(ybus(k,nn))*abs(V(nn))*cos(-angle(ybus(k,nn))-th(nn)+th(k));
end
end
end
end
end
end
% Calculate J22
for k=2:nbus-pvbus
for n=2:nbus-pvbus
if busdata(n,2)== 3 && busdata(k,2)== 3
if k~=n
J22(k-1,n-1)= abs(V(k))*abs(ybus(k,n))*sin(-angle(ybus(k,n))-th(n)+th(k));
end
if k==n
J22(k-1,k-1) = -(J11(k-1,k-1))./abs(V(k-1))-(2*abs((V(k))))*imag(ybus(k,k));
end
end
end
end
% Update J
display(['*********************',' Iteration # ',num2str(iteration),'***********************'])
display('********************* Jacobian Matrix *********************')
J = [J11 J12;J21 J22]
delta_theta_V = inv(J)*deltaPQ;
theta_V = theta_V + delta_theta_V;
% Update V and Theata
ite = 1;
for n=2:nbus
if busdata(n,2)==3
V(n) = theta_V(nbus-1+ite);
ite = 1+ite;
end
end
th = [0;theta_V(1:nbus-1)];
display('********************* Voltage Angle and Voltage Magnitude *********************')
Vmag = V % Bus Voltages
Ang = 180/pi*th % Bus Voltage Angles in Degree
Voltage = V.*cos(th)+i.*V.*sin(th) % Bus Voltages in phasor form
Va = V(2:end);
tha = th(2:end);
i=0;
for k=1:nbus
if busdata(k,2)== 2
Va(k-1-i) = [];
i=i+1;
end
end
theta_V = [tha;Va];
toler = max(abs(V-Vpe));
if toler<0.00001
for n=1:nbus
if type(n)==1 | 2
QQ=0;
for i=1:nbus;
QQ = QQ + abs(V(n))*abs(V(i))*abs(ybus(n,i))*sin(th(n)-th(i)-angle(ybus(n,i)));
end
PP=0;
for i=1:nbus
PP = PP + abs(V(n))*abs(V(i))*abs(ybus(n,i))*cos(th(n)-th(i)-angle(ybus(n,i)));
end
P(n)=PP;
Q(n)=QQ;
end
end
display('********************** Final Results **********************')
iteration % Total iterations
VoltagePhasor = V.*cos(th)+i.*V.*sin(th) % Final bus Voltages in phasor form
VoltageMagnitude = V % Final Bus Voltage Magnitude (per unit)
VoltageAngle_Degree = 180/pi*th % Final Bus Voltage Angles in Degree
display('*************************** END ***************************')
return
end
Vpe = V;
end
%%
douglas-2013
Newton-RaphsonMethod1.docx
0
