Homework SourseGenerators G1 100 MVA 138 kV G2 200 MVA 150 kV Transformers
Generators
G1: 100 MVA, 13.8 kV
G2: 200 MVA, 15.0 kV
Transformers
T1: 100 MVA, 13.8 kV/230 kV, XT1=0.1 pu
T2: 200 MVA, 15.0 kV/230 kV, XT2=0.1 pu
Transmission Lines
All lines 230 kV, z = 0.08+j0.5 ohm/km
Line lengths are
Line-1=15 km,
Line-2=50 km,
Line-3=40 km,
Line-4=15 km,
Line-5=50 km
System Base Quantities
Sbase = 100 MVA (Three-phase)
Vbase = 13.8 Kv (Line-to-line) in the zone of G1
Generator, line, transformer data information and one-line diagram of 7 bus system are given above.
Besides, all transformers have same turn ratio which is 1. Using the given system, solve the system for power flow analysis with Newton Raphson Method.
MATPOWER simulation package will be used to solve power flow problem for the given system. You can download the MATPOWER from the following link.
http://www.pserc.cornell.edu/matpower/
Please include data file and the power flow results.
If you use any calculations during the preparation of the data file, please also include them in the report.
G1 G2 180 MW V1 13,8kV BUS 1 BUS 7 IV71-15.0kV T2 T1 BUS 2 BUS 5 Line-5 30 MVAr 30 MVAr 50 MW Line-1 50 MW Line-4 BUS 3 Line-3 BUS 4 Line-2 30 MVAr 50 MW 50 MW 30 MVAr 50 MW BUS 30 MVAr 6 The on-line diagram of 7 bus systemSolution
data in the type
% b itr rtr nrtr trl sz
% ------------------------------------------------------------------------------
5 110 0 0 7 0 0 0 0 0
% pq g tol
% ------------------------------------------------------------------------------
3 2 .000001 0 0 0 0 0 0 0
%reg.trsf.data
%From To R X B/2
% ------------------------------------------------------------------------------
%1 3 .00 .00 0 0 0
%n-reg. trsf.data
%From To R X B/2
% ------------------------------------------------------------------------------
%2 5 .000625 .003125 0 0 0
%transmission line data
%From To R X B/2
% ------------------------------------------------------------------------------
1 2 .0200 .0600 0 0.030 1 0 0 0
1 3 .0800 .2400 0 0.025 1 0 0 0
2 3 .0600 .1800 0 0.020 1 0 0 0
2 4 .0600 .1800 0 0.020 1 0 0 0
2 5 .0400 .1200 0 0.015 1 0 0 0
3 4 .0000 .1000 0 0.000 1 0 0 0
4 5 .0800 .2400 0 0.025 1 0 0 0
%shunt impedance data
% Bus R X
% ------------------------------------------------------------------------------
%3 0 .1 0 0 0 0
% P - Q Load Data
% Bus MW MVAR
% ------------------------------------------------------------------------------
3 0.45 0.15 0.00 .00 0 0 0 0 0
4 0.40 0.05 0.00 .00 0 0 0 0 0
5 0.60 0.10 0.00 .00 0 0 0 0 0
% Generation Data
% Bus MW Qgmx Qgmn Vsp Vmax Vmin
% ------------------------------------------------------------------------------
1 1.2 .5 -3.0 1.0600 1.150 1.000 0 0 0
2 .20 .20 -3.0 1.0000 1.150 1.100 0 0 0
program:
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 pqdata1(i,1)==pvdata(i,1),
loaddata(i,10)=pqdata1(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;
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))));
itr=itr+1;
if tol<=0.0001,
disp(\'the system is converging in the no.of iterations= \');
disp(itr);
end
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);
