Homework SourseFor the given system data construct Ybus matrix Solve for un
Solution
clc
clear all;
load data.m;
bus=data(1,:);
nb=bus(1);
nl=bus(3)+bus(5);
tibus=data(2,:);
Npv=tibus(2);
npq=tibus(1);
Npq=nb-Npv;
linedata=data(3:nl+2,1:end);
sn=bus(6);
if sn~=0,
bu=data(nl+3:nl+2+sn,1:end);
else
bu=zeros(1,end);
end
pqdata1=zeros(nb,1);
pqdata=data(nl+3+sn:(nl+npq++sn+2),(1:end));
pvdata=data((nl+npq+sn+3):(nl+npq+sn+Npv+2),(1:end));
for i=1:10,
pqdata1(pqdata(:,1),i)=pqdata(:,i);
end
for j=1:nb,
if pqdata1(j,1)==0,
pqdata1(j,1)=j;
end
end
pqdata1(:,5)=1;
pqdata1(:,6)=max(pvdata(:,6));
pqdata1(:,7)=min(pvdata(:,7));
pqdata1(:,8)=pqdata1(:,3);
pqdata(:,3)=0;
for i=1:npq,
if i<=Npv,
loaddata=pvdata;
end
end
for i=1:Npv,
if pqdata(i,1)==pvdata(i,1),
loaddata(i,10)=pqdata(i,2);
end
end
loaddata=[loaddata; pqdata1(Npv+1:nb,1:end)];
disp(\'---------------------------------------------% given data%-------------------------------------------------------------------------------\');
disp(\'----------------------------------------------%bus data%---------------------------------------------------------------------------------\');
disp(\' b itr rtr nrtr trl sz SVC Pd\');
disp(bus);
disp(\'-----------------------------------------------%transmission bus data%--------------------------------------------------------------------\');
disp(\' pq g tol\');
disp(tibus);
disp(\'------------------------------------------------%line data%--------------------------------------------------------------------------------\');
disp(\' From To R X G/2 B/2 TAP limits tmin tmax\');
disp(linedata);
disp(\'-----------------------------------------------%shunt admitance data%-----------------------------------------------------------------------\');
disp(\'Bus R X\')
disp(bu);
disp(\'------------------------------------------------%loadf data%--------------------------------------------------------------------------------\');
disp(\' Bus MW Qgmx Qgmn Vsp Vmax Vmin\');
disp(loaddata);
itr=0;
a=linedata(:,7); %---------Ybus caluculation------------%
r=linedata(:,3);
x=linedata(:,4);
z=r+1i*x;
b=zeros(nb,1);
if bu(:,1)~=0,
b(bu(:,1))=1./bu(:,4);
end
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);
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);
if max(loaddata(:,2))>50,
base=100;
else
base=1;
end
qlmx=loaddata(:,3)/base;
qlmn=loaddata(:,4)/base;
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))/base;
end
for i=Npv+1:nb,
Qspec(i-(Npv),1)=loaddata(i,8)/base;
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 in per unit..........,\');
disp(pl);
if base==100,
disp(\'power loss in MWs..................,\');
pl=pl*base;
disp(pl);
end
sf=zeros(nl,3);
for i=1:nl,
sf(i,1)=fb(i);
sf(i,2)=tb(i);
sf(i,3)=vrect(fb(i))*((conj(vrect(fb(i)))-conj(vrect(tb(i))))*a(i)*(-conj(Y(fb(i),tb(i))))+conj(vrect(fb(i)))*conj(complex(0,0.5*linedata(i,6))));
end
sf=conj(sf);
disp(\'------------------------%forward power flow%-------------------------------------------\');
disp(\' from to power flow\');
disp(sf);
sr=zeros(nl,3);
for i=1:nl,
sr(i,1)=tb(i);
sr(i,2)=fb(i);
sr(i,3)=vrect(tb(i))*((conj(vrect(tb(i)))-conj(vrect(fb(i))))*a(i)*(-conj(Y(tb(i),fb(i))))+conj(vrect(tb(i)))*conj(complex(0,0.5*linedata(i,6))));
end
sr=conj(sr);
disp(\'------------------------%reverse power flow%-------------------------------------------\');
disp(\' from to power flow\');
disp(sr);
