Please solve all parts Consider the power system shown below. All data are in pe
ID: 2291538 • Letter: P
Question
Please solve all parts
Consider the power system shown below. All data are in per-unit. Circuit parameters are per unit admittances. Generators are located at buses 1, 2 and 5 1-j10 -j20 1-j12 1-j10 Additional data P,2-P,,-2.0 S3 S42.5+j0.0 1. Find the bus admittance matrix Ybus for this system. Find the matrices G and B GS G3 where Ybus-G+jB 2. Identify each bus type, i.e. slack, PV or PQ (justify your choice) 3. What are the unknown load flow state variables 4. Write the load flow equations for this system in terms of the unknown load flow state variables and circuit parameters 5. Use the Newton-Raphson algorithm to find the load flow state variables. For a 6. 7. 8. 9. stopping criterion use ?=10. You can use MATLAB at least for matrix operations (invert J, find mismatch, etc) or you can write a short programm. Compute the real and reactive power generated at the slack bus Compute the reactive power generated at the PV buses Compute all line flows and shunt flows Summarize results on a single-line diagram (if it is crowded, do one for real power and a second one for reactive power)Explanation / Answer
Answer:
MATLAB CODE FOR NEWTON RAPHSON LOAD FLOW METHOD FOR 5BUS SYSTEM:
clc
clear all
load linedata5c.m
A=linedata5c;
nb=A(1,1);
t=A(1,2);
rtr=A(1,3);
nrtr=A(1,4);
ntr=A(1,5);
sz=A(1,6);
pq=A(2,1);
g=A(2,2);
nl=rtr+ntr+nrtr;
m=nl+2+sz;
Ld=A(nl+sz+3:m+pq,:);
Gd=A(nl+sz+pq+3:m+pq+g,:);
for i=1:nl
sb(i)=A(i+2,1);
rb(i)=A(i+2,2);
R(i)=A(i+2,3);
X(i)=A(i+2,4);
B(i)=A(i+2,6)/2;
bx(i)=complex(0,B(i));
TP(i)=A(i+2,7);
Z(i)=complex(R(i),X(i));
y(i)=1/Z(i);
end
ybus=zeros(nb,nb);
for j=1:nl
ybus(sb(j),rb(j))=ybus(sb(j),rb(j))-y(j)/TP(j);
ybus(rb(j),sb(j))=ybus(sb(j),rb(j));
end
ybus;
for k=1:nb
for l=1:nl
if sb(l)==k
ybus(k,k)=ybus(k,k)+y(l)/(TP(l)^2)+bx(l);
else if rb(l)==k
ybus(k,k)=ybus(k,k)+y(l)+bx(l);
end
end
end
end
ybus;
Y1=zeros(nb,nb);
s=A(nl+3:m,:);
sh=s(:,1);
for i=1:sz
if ((sz~=0)&(s(i,1)==sh(i,1)))
Y1(sh(i,1),sh(i,1))=Y1(sh(i,1),sh(i,1))+1/complex(s(i,3),s(i,4));
end
end
Y1;
if sz~=0
y1=[ybus+Y1];
else
y1=ybus;
end
y1;
G=real(y1);
B1=imag(y1);
Y=abs(ybus);
YA=angle(ybus);
VM=ones(nb,1);
for k=1:g
% if Gd(k,1)==gd(k,1)
VM(k)=A(k+m+pq,5);
% end
end
VM;
DL=zeros(nb,1);
iter=0
while iter<=t
Pcal=zeros(nb,1);
Qcal=zeros(nb,1);
for i=1:nb
for k=1:nb
Pcal(i) = Pcal(i) + VM(i)* VM(k)*(G(i,k)*cos(DL(i)-DL(k)) + B1(i,k)*sin(DL(i)-DL(k)));
Qcal(i) = Qcal(i) + VM(i)* VM(k)*(G(i,k)*sin(DL(i)-DL(k)) - B1(i,k)*cos(DL(i)-DL(k)));
end
end
Pcal;
Qcal;
%Mismatch Calculation
Ps=zeros(nb,1);
bl=Ld(:,1);
for i=1:g
Ps(i)=Ps(i)+Gd(i,2);
end
for i1=1:pq
Ps(bl(i1))=Ps(bl(i1))-Ld(i1,2);
end
Qs=zeros(nb,1);
for i1=1:pq
Qs(bl(i1))=Qs(bl(i1))-Ld(i1,3);
end
Ps;
Qs;
for q=1:nb
DP(q)=Ps(q)-Pcal(q);
DQ(q)=Qs(q)-Qcal(q);
end
DP;
DQ;
DPA=DP(2:nb);
DQA=DQ(g+1:nb);
DPQ=[DPA';DQA'];
if max(abs(DPQ))<=0.00001
disp('----------------converged------------------')
Pcal
Qcal
VM
DL
Ploss=sum(Pcal)
Qloss=sum(Qcal)
Ygb2=zeros(nb,nb);
for xk=1:nl
Ygb2(sb(xk),rb(xk))=bx(xk);
Ygb2(rb(xk),sb(xk))=Ygb2(sb(xk),rb(xk));
end
Ygb2;
g2=real(Ygb2);
b2=imag(Ygb2);
for pk=1:nb
for qk=1:nb
if pk~=qk
PF(pk,qk)=-VM(pk)*VM(pk)*Y(pk,qk)*cos(YA(pk,qk))+VM(pk)*VM(qk)*Y(pk,qk)*cos(DL(pk)-DL(qk)-YA(pk,qk));
QF(pk,qk)=VM(pk)*VM(pk)*Y(pk,qk)*sin(YA(pk,qk))+VM(pk)*VM(qk)*Y(pk,qk)*sin(DL(pk)-DL(qk)-YA(pk,qk))-VM(pk)*VM(pk)*b2(pk,qk);
end
end
end
PF %required power flow
QF
break
else
% Formation of Jacobian Matrix
for m1=2:nb
for n1=2:nb
if m1==n1
DPDD(m1,m1)=-Qcal(m1)-B1(m1,m1)*VM(m1)*VM(m1);
else
DPDD(m1,n1)=VM(m1)*VM(n1)*(G(m1,n1)*sin(DL(m1)-DL(n1))-B1(m1,n1)*cos(DL(m1)-DL(n1)));
end
end
end
for m2=2:nb
for n2=g+1:nb
if m2==n2
DPDVM(m2,m2)=Pcal(m2)/VM(m2)+G(m2,m2)*VM(m2);
else
DPDVM(m2,n2)=VM(m2)*(G(m2,n2)*cos(DL(m2)-DL(n2))+B1(m2,n2)*sin(DL(m2)-DL(n2)));
end
end
end
for m3=g+1:nb
for n3=2:nb
if m3==n3
DQDD(m3,m3)=Pcal(m3)-G(m3,m3)*VM(m3)*VM(m3);
else
DQDD(m3,n3)=-VM(m3)*VM(n3)*(G(m3,n3)*cos(DL(m3)-DL(n3))+B1(m3,n3)*sin(DL(m3)-DL(n3)));
end
end
end
for m4=g+1:nb
for n4=g+1:nb
if m4==n4
DQDVM(m4,m4)=Qcal(m4)/VM(m4)-B1(m4,m4)*VM(m4);
else
DQDVM(m4,n4)=VM(m4)*(G(m4,n4)*sin(DL(m4)-DL(n4))-B1(m4,n4)*cos(DL(m4)-DL(n4)));
end
end
end
J1=DPDD(2:nb,2:nb);
J2=DPDVM(2:nb,g+1:nb);
J3=DQDD(g+1:nb,2:nb);
J4=DQDVM(g+1:nb,g+1:nb);
J=[J1 J2;
J3 J4];
JI=inv(J);
DDDV=JI*DPQ;
%DDL=zeros(nb,1);
%DVM=zeros(nb,1);
DDL=DDDV(1:nb-1)
DVM=DDDV(nb:2*(nb)-1-g)
for a=1:nb-1
DL(a+1)=DL(a+1)+DDL(a);
end
for b=1:nb-g
VM(b+g)=VM(b+g)+DVM(b);
end
DL;
VM;
end
iter=iter+1
end
PQ=complex(PF,QF)
% Bus Power Injections..
Si=zeros(nl,1)
for i=1:nl
for l=1:nb
for k1=1:nb
if l~=k1
Si(i)=Si(i)+PQ(l,k1);
end
end
end
end
Si
Pi = real(Si);
Qi = -imag(Si);
for i=1:nl
SB=sb(i);
RB=rb(i);
end
disp('***********************************************');
disp(' SB RB Pi Qi ');
disp('**************************************************');