Homework SourseKy 020 81 030 Solutionclc clear all load busdat load lineda
Solution
clc;
clear all;
load bus.dat;
load linedata.asv;
load bu.dat;
itr=0;
a=linedata(:,7); %---------Ybus caluculation------------%
nb=bus(1,1);
nl=bus(1,3)+bus(1,5);
r=linedata(:,3);
x=linedata(:,4);
z=r+1i*x;
b=zeros(nb,1);
c=1i.*b;
disp(\'b\');
disp(b)
y=1./z;
disp(\'y\');
disp(y)
fb=linedata(:,1);
tb=linedata(:,2);
v=ones(nb,1);
d=zeros(nb,1);
Y=zeros(nb,nb);
for i=1:nl
Y(fb(i),tb(i))=Y(fb(i),tb(i))-(y(i)/conj(a(i)));
Y(tb(i),fb(i))=Y(fb(i),tb(i));
end
for i=1:nb
for j=1:nl
if fb(j)==i
Y(i,i)=Y(i,i)+y(j)/((a(j))^2)+complex(0,linedata(j,6));
else if tb(j)==i
Y(i,i)=Y(i,i)+y(j)+complex(0,linedata(j,6));
end
end
end
Y(i,i)=Y(i,i)-c(i);
end
disp(\'Y\');
disp(Y)
G=real(Y);
B=imag(Y);
load tibus.dat;
load loaddata.asv;
Npq=tibus(1,1);
Npv=tibus(1,2);
g=2*Npq+Npv;
bus=ones(nb,1);
for i=Npv+1:nb,
bus(i,1)=2;
end
for i=1:nb
v(i,1)=loaddata(i,5);
d(i,1)=0;
end
disp(\'v\');
disp(v)
vmx=loaddata(:,6);
vmn=loaddata(:,7);
qlmx=loaddata(:,3)/100;
qlmn=loaddata(:,4)/100;
tol=10;
%--------NEWTON RAPHSON METHOD----------
while (tol>0.0001)
vprev=v;
P=zeros(nb,1);
Q=zeros(nb,1);
for i=1:nb
for k=1:nb
P(i,1)=P(i,1)+(v(i,1))*(G(i,k)*cos(d(i,1)-d(k,1))+B(i,k)*sin(d(i,1)-d(k,1)))*v(k,1);
Q(i,1)=Q(i,1)+(v(i,1))*(G(i,k)*sin(d(i,1)-d(k,1))-B(i,k)*cos(d(i,1)-d(k,1)))*v(k,1);
end
if bus(i)==1,
if (Q(i)>qlmx(i) || Q(i)<qlmn(i)),
if(Q(i)<qlmn(i)),
Q(i)=qlmn(i);
else
Q(i)=qlmx(i);
end
bus(i)=2;
else
bus(i)=1;
end
end
end
disp(\'bus\')
disp(bus)
disp(\'P\')
disp(P)
disp(\'Q\')
disp(Q)
Pspec=zeros(nb-1,1);
Qspec=zeros(Npq,1);
for i=2:nb
Pspec(i-1,1)=(loaddata(i,2)-loaddata(i,10))/100;
end
for i=Npv+1:nb,
Qspec(i-(Npv),1)=loaddata(i,8)/100;
end
delP=zeros(nb-1,1); %--------mismatch calculation----------------%
delQ=zeros(Npq,1);
for i=2:Npv
delP(i-1,1)=Pspec(i-1,1)-P(i,1);
end
for i=Npv+1:nb
delP(i-1,1)=(-1)*Pspec(i-1,1)-P(i,1);
delQ(i-(Npv),1)=(-1)*Qspec(i-(Npv),1)-Q(i,1);
end
disp(\'delP\');
disp(delP)
disp(\'delQ\');
disp(delQ)
del=zeros(g-1,1);
for i=1:g-1
if (i<nb)
del(i,1)=delP(i,1);
end
if (i>=nb)
del(i,1)=delQ(i-(nb-1),1);
end
end
disp(\'del\');
disp(del)
J11=zeros(nb-1,nb-1); %-----------jacobian matrix calculation------------%
J12=zeros(nb-1,Npq);
J21=zeros(Npq,nb-1);
J22=zeros(Npq,Npq);
for i=2:nb
for k=2:nb
if(k==i)
for l=1:nb
if(l~=i)
J11(i-1,k-1)=J11(i-1,k-1)+v(i,1)*v(l,1)*(-G(i,l)*sin(d(i,1)-d(l,1))+B(i,l)*cos(d(i,1)-d(l,1)));
end
end
else
J11(i-1,k-1)=(v(i,1)*v(k,1)*(G(i,k)*sin(d(i,1)-d(k,1))-B(i,k)*cos(d(i,1)-d(k,1))));
end
end
end
for i=2:nb
for k=Npv+1:nb
if(i==k)
for l=1:nb
if(l~=i)
J12(i-1,k-(Npv))=J12(i-1,k-(Npv))+v(l,1)*(G(i,l)*cos(d(i,1)-d(l,1))+B(i,l)*sin(d(i,1)-d(l,1)));
end
end
J12(i-1,k-(Npv))=J12(i-1,k-(Npv))+2*v(i,1)*G(i,i);
else
J12(i-1,k-(Npv))=v(i,1)*(G(i,k)*cos(d(i,1)-d(k,1))+B(i,k)*sin(d(i,1)-d(k,1)));
end
end
end
for i=Npv+1:nb
for k=2:nb
if(i==k)
for l=1:nb
if(l~=i)
J21(i-(Npv),k-1)=J21(i-(Npv),k-1)+v(i,1)*v(l,1)*(G(i,l)*cos(d(i,1)-d(l,1))+B(i,l)*sin(d(i,1)-d(l,1)));
end
end
else
J21(i-(Npv),k-1)=(v(i,1)*v(k,1)*(-G(i,k)*cos(d(i,1)-d(k,1))-B(i,k)*sin(d(i,1)-d(k,1))));
end
end
end
for i=Npv+1:nb
for k=Npv+1:nb
if(k==i)
for l=1:nb
if(l~=i)
J22(i-(Npv),k-(Npv))=J22(i-(Npv),k-(Npv))+(v(l,1)*(G(i,l)*sin(d(i,1)-d(l,1))-B(i,l)*cos(d(i,1)-d(l,1))));
end
end
J22(i-(Npv),k-(Npv))=J22(i-(Npv),k-(Npv))-2*v(i,1)*B(i,i);
else
J22(i-(Npv),k-(Npv))=v(i,1)*(G(i,k)*sin(d(i,1)-d(k,1))-B(i,k)*cos(d(i,1)-d(k,1)));
end
end
end
disp(\'J11\');
disp(J11)
disp(\'J12\');
disp(J12)
disp(\'J21\');
disp(J21)
disp(\'J22\');
disp(J22)
J=[J11 J12;J21 J22];
disp(\'J\');
disp(J)
Ji = inv(J);
Vcn=(inv(J))*del; %---------change in load angles and voltages----------%
disp(\'Vcn\');
disp(Vcn)
deld=zeros(nb-1,1);
delv=zeros(Npq,1);
for i=1:nb-1
deld(i,1)=Vcn(i,1);
end
disp(\'deld\');
disp(deld)
for i=nb:g-1
delv(i-(nb-1),1)=Vcn(i,1);
end
disp(\'delv\');
disp(delv)
dnew=zeros(nb-1,1);
vnew=zeros(Npq,1);
for i=2:nb %-----------updating the voltages and load angles-----------%
d(i,1)=d(i,1)+deld(i-1,1);
end
disp(\'d\');
disp(d)
for i=Npv+1:nb
v(i,1)=v(i,1)+delv(i-Npv,1);
end
disp(\'v\');
disp(v)
tol=max(abs(abs(v(:,1))-abs(vprev(:,1))));
if tol<=0.0001,
disp(\'the system is converging in the no.of iterations= \');
disp(itr);
end
itr=itr+1;
disp(\'tol\');
disp(tol)
end
p1=0;q1=0;
for k=1:nb
P(1,1)=p1+(v(1,1))*(G(1,k)*cos(d(1,1)-d(k,1))+B(1,k)*sin(d(1,1)-d(k,1)))*v(k,1);
Q(1,1)=q1+(v(1,1))*(G(1,k)*sin(d(1,1)-d(k,1))-B(1,k)*cos(d(1,1)-d(k,1)))*v(k,1);
p1=P(1);
q1=Q(1);
end
disp(\'P\');
disp(P)
disp(\'Q\');
disp(Q)
d=d*(180/pi);
disp(\'VM\');
disp(v);
disp(\'delta_bus load angle\');
disp(d);
vrect=zeros(nb,1);
for i=1:nb,
vrect(i)=v(i)*(cosd(d(i))+1i*sind(d(i)));
end
pls=zeros(11,4);
plss=zeros(11,4);
I0=zeros(nl,1);
Ic=zeros(nl,1);
sflow=zeros(nl,1);
I=zeros(nl,1);
ploss=zeros(nl,1);
for i=1:nl,
for j=1:nb,
for k=1:nb,
if fb(i)==j && tb(i)==k,
I0(i)=1i*b(j)*vrect(j);
Ic(i)=v(j)*complex(0,linedata(i,6));
I(i)=(conj(vrect(j))-conj(vrect(k)))*(-(conj(Y(j,k))));
ploss(i,1)=(abs(I(i))^2)*real(1./(-(conj(Y(j,k)))));
pls(i,1)=i;
pls(i,2)=j;
pls(i,3)=k;
pls(i,4)=ploss(i,1);
end
end
end
end
disp(\'----------------%power loss%-------------------\');
disp(\'tr.line no. bus(from) bus(to) powerloss\');
disp(pls)
pl=0;
for i=1:nl,
pl=pl+pls(i,4);
end
disp(\'total powerloss\');
disp(pl);
disp(\'power loss in MWs........\');
pl=pl*100;
disp(pl);
