y b u s

code for y-bus matrix

[KALYAN ACHARJYA]

Please complete the question?

Answer 1

[Vyshnavi]

#exp-1 calculation of inductance and capacitance

type=input('enter the type of power system')

disp('1:1 phase');

disp('2:3 phase double circuit');

disp('3:3 phase symmetrically placed');

disp('4:3 phase unsymmetrically placed');

if type==1

d=input('enter distance b/w 2 wires');

r=input('enter the radius of the wire');

r1=0.7788*r;

L=0.2.*log(d./r1);

fprintf('inductance is %f mH/km\n',L);

C=0.0555./log(d./r1);

fprintf('capacitance is %f uF/km\n',C);

elseif type==2

l=input('enter the distance b/w 2 transposed line in m');

d=input('Distance two conductors in line in m');

r=input('enter the radius of the wire');

r1=0.7788*r;

m=power((power(l,2)+power(d,2)),0.5);

n=power((power(2.*d,2)+power(l,2)),0.5);

x=power(2,1/6).*power((d./r1),1/2).*power(m/n,1/3);

y=power(2,1/3).*(d./r).*power(m/n,2/3);

L=0.2.*log(x);

C=0.1111.*(1./log(y));

fprintf('inductance is %f mH/km\n',L);

fprintf('capacitance is %f uF/km\n',C);

elseif type==3

d=input('enter distance b/w 2 wires');

r=input('enter the radius of the wire');

r1=0.7788*r;

L=0.2.*log(d./r1);

C=0.0555./log(d./r);

fprintf('inductance is %f mH/km\n',L);

fprintf('capacitance is %f uF/km\n',C);

elseif type==4

d12=input('enter distance b/w 1 & 2 wires in m');

d23=input('enter distance b/w 2 & 3 wires in m');

d31=input('enter distance b/w 1 & 3 wires in m');

r=input('enter the radius of the wire');

r1=0.7788*r;

L=0.2.*log(d./r1);

C=0.0555./log(d./r);

fprintf('inductance is %f mH/km\n',L);

fprintf('capacitance is %f uF/km\n',C);

end

#exp-2 bus admittance matrix

n=input('enter the number of buses ');

for q=1:n

for r=q+1:n

fprintf('enter the impedance value between %d-%d',q,r);

Z(q,r)=input('=');

if (Z(q,r)==0)

y(q,r)=0;

else

y(q,r)=inv(Z(q,r));

end

y(r,q)=y(q,r);

fprintf('enter the half line charging admittance');

X(q,r)=input('=');

X(r,q)=X(q,r);

end

end

tr=zeros(n,n);

for s=1:n

fprintf('enter the self admittance of the bus %d',s);

u(s)=input('=');

end

ybus=zeros(n,n);

for a=1:n

for b=1:n

if (a==b)

for c=1:n

ybus(a,a)=ybus(a,a)+y(a,c)+X(a,c);

end

else

ybus(a,b)=-y(b,a);

end

end

ybus(a,a)=ybus(a,a)+u(a);

end

for r=1:n

for h=1:n

if (r==h)

ybus(r,r)=ybus(r,r)+tr(r,r);

else

ybus(r,h)=-(y(r,h)+tr(r,h));

end

end

end

ybus

#exp-3 bus impedance matrix

disp('Formation of Impedence Matrix using Bus building Algorithm')

n = input('Enter total number of buses including reference bus=');

zbus = zeros(n,n);

t=1;

while t==1

zbus;

S = menu ('Specify case no', 'New bus to reference bus','Existing bus to new bus','Between existing buses','Existing buses to reference bus','Print','Quit');

switch S

case{1}

zb = input('Enter impedence value=');

zbus = zb;

case{2}

k = input('Enter the starting bus number=');

n = input('Enter new bus number=');

zb = input('Enter impedence value=');

for i= 1:n

if i==n

zbus(n,n) = zbus(k,k)+zb;

else

zbus(i,n) = zbus(i,k);

zbus(n,i) = zbus(k,i);

end

end

case{3}

a= input('Enter first bus number=');

b = input('Enter second bus number=');

zb =input('Enter impedence value=');

m1= zb+zbus(a,a)+zbus(b,b)-(2*zbus(a*b));

ztemp = (1/m1)*((zbus(:,a))-(zbus(:,b)))*(((zbus(a,:))-zbus(b,:)));

zbus = zbus-ztemp;

case{4}

k = input('Enter the old bus no=');

zb = input('Enter impedence value=');

m2 = zbus(k,k)+zb;

ztemp = (1/m2)*zbus(:,k)*zbus(k,:);

zbus = zbus-ztemp;

case{5}

zbus

case{6}

'end program';

t=0;

end

end

#exp-4 load flow analysis using gauss seidel method

%mainfunction file

busdata=[1 1 1.025 0 0 0 0 0

2 0 1.0 0 400 200 0 0

3 1 1.0 0 0 0 300 0 ];

linedata=[1 2 0 0.025 0 1

1 3 0 0.05 0 1

2 3 0 0.025 0 1];

%Lfybus file:

% Y Bus admittance matrix

j=sqrt(-1); i = sqrt(-1);

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); % Shunt Admittance, B/2

Z = R + j*X; % Z matrix

y= 1./Z; %branch admittance

nbranch=length(linedata(:,1)); % no. of branches

nbus = max(max(fb), max(tb)); % no. of buses

% Forming the Y Bus Matrix

for n = 1:nbranch

Ybus=zeros(nbus,nbus); % initialize Ybus to zero

% Formation of the off diagonal elements

for k=1:nbranch

Ybus(fb(k),tb(k))=Ybus(fb(k),tb(k))-y(k);

Ybus(tb(k),fb(k))=Ybus(fb(k),tb(k));

end

end

% Formation of the diagonal elements

for n=1:nbus

for k=1:nbranch

if fb(k)==n || tb(k)==n

Ybus(n,n) = Ybus(n,n)+y(k) + b(k);

else, end

end

end

%Lfgauss file:

% Gauss-Seidel method

basemva = 100; %Base MVA

tolerance = 0.001; %Tolerance

mi = 100; %Maximum Iterations

af = 1.8; %Acceleration factor

% Keys for check purposes

Vm=0; delta=0; yload=0; deltad =0;

nbus = length(busdata(:,1));

for k=1:nbus

n=busdata(k,1);

bt(n)=busdata(k,2);

Vm(n)=busdata(k,3);

delta(n)=busdata(k, 4);

Pd(n)=busdata(k,5);

Qd(n)=busdata(k,6);

Pg(n)=busdata(k,7);

Qg(n) = busdata(k,8);

if Vm(n) <= 0

Vm(n) = 1.0; V(n) = 1 + j*0;

else delta(n) = pi/180*delta(n);

V(n) = Vm(n)*(cos(delta(n)) + j*sin(delta(n)));

P(n)=(Pg(n)-Pd(n))/basemva;

Q(n)=(Qg(n)-Qd(n))/basemva;

S(n) = P(n) + j*Q(n);

end

DV(n)=0;

end

num = 0; AcurBus = 0; converge = 1;

Vc = zeros(nbus,1)+j*zeros(nbus,1); Sc = zeros(nbus,1)+j*zeros(nbus,1);

iter=0;

maxerror=10;

sumc1=0;

sumc2=0;

sumc3=0;

sumc4=0;

while maxerror >= tolerance && iter <= mi

iter=iter+1;

for n = 1:nbus

YV = 0+j*0;

for L = 1:nbranch

if fb(L) == n, k=tb(L);

YV = YV + Ybus(n,k)*V(k);

elseif tb(L) == n, k=fb(L);

YV = YV + Ybus(n,k)*V(k);

end

end

Sc = conj(V(n))*(Ybus(n,n)*V(n) + YV) ;

Sc = conj(Sc);

DP(n) = P(n) - real(Sc);

DQ(n) = Q(n) - imag(Sc);

if bt(n) == 1

S(n) =Sc; P(n) = real(Sc); Q(n) = imag(Sc); DP(n) =0; DQ(n)=0;

Vc(n) = V(n);

elseif bt(n) == 2

Q(n) = imag(Sc); S(n) = P(n) + j*Q(n);

Qgc = Q(n)*basemva + Qd(n) ;

end

if bt(n) ~= 1

Vc(n) = (conj(S(n))/conj(V(n)) - YV )/ Ybus(n,n);

else, end

if bt(n) == 0

V(n) = V(n) + af*(Vc(n)-V(n));

elseif bt(n) == 2

VcI = imag(Vc(n));

VcR = sqrt(Vm(n)^2 - VcI^2);

Vc(n) = VcR + j*VcI;

V(n) = V(n) + af*(Vc(n) -V(n));

end

end

maxerror=max( max(abs(real(DP))), max(abs(imag(DQ))) );

if iter == mi && maxerror > tolerance

fprintf('\nWARNING: Iterative solution did not converge after ')

fprintf('%g', iter), fprintf(' iterations.\n\n')

fprintf('Press Enter to terminate the iterations and print the results \n')

converge = 0; pause, else, end

end

if converge ~= 1

tech= (' ITERATIVE SOLUTION DID NOT CONVERGE');

else

tech=('POWER FLOW SOLUTION: GAUSS SEIDAL METHOD');

end

k=0;

for n = 1:nbus

Vm(n) = abs(V(n)); deltad(n) = angle(V(n))*180/pi;

if bt(n) == 1

S(n)=P(n)+j*Q(n);

Pg(n) = P(n)*basemva + Pd(n);

Qg(n) = Q(n)*basemva + Qd(n) ;

k=k+1;

Pgg(k)=Pg(n);

elseif bt(n) ==2

k=k+1;

Pgg(k)=Pg(n);

S(n)=P(n)+j*Q(n);

Qg(n) = Q(n)*basemva + Qd(n) ;

end

sumc1 = sumc1 + Pg(n);

sumc2 = sumc2+ Qg(n);

sumc3 = sumc3 + Pd(n);

sumc4 = sumc4 + Qd(n);

yload(n) = (Pd(n)- j*Qd(n))/(basemva*Vm(n)^2);

end

busdata(:,3)=Vm'; busdata(:,4)=deltad';

%Busout file:

disp(tech)

fprintf(' %g Iterations \n\n', iter)

head =[' BusNo Voltage(Mag) Angle(degree)----Load (MW,MVar)-----Generation(MW,MVar)'];

disp(head)

for n=1:nbus

fprintf(' |%5g', n), fprintf(' |%7.3f', Vm(n)),

fprintf(' |%8.3f', deltad(n)), fprintf(' |%9.3f', Pd(n)),

fprintf(' |%9.3f', Qd(n)), fprintf(' |%9.3f', Pg(n)),

fprintf(' |%9.3f \n', Qg(n)),

end

fprintf(' \n'), fprintf(' Total ')

fprintf(' %9.3f', sumc3), fprintf(' %9.3f', sumc4),

fprintf(' %9.3f', sumc1), fprintf(' %9.3f\n\n',sumc2),

%Lineflow file:

SLT = 0;

fprintf('\n')

fprintf(' Line Flow and Losses \n\n')

fprintf(' --Line-- Power at bus & line flow --Line loss-- \n')

fprintf(' from to MW Mvar MVA MW Mvar \n')

for n = 1:nbus

busprt = 0;

for L = 1:nbranch

if busprt == 0

fprintf(' \n'), fprintf('%6g', n), fprintf(' %9.3f', P(n)*basemva)

fprintf('%9.3f', Q(n)*basemva), fprintf('%9.3f\n', abs(S(n)*basemva))

busprt = 1;

else, end

if fb(L)==n k = tb(L);

In = (V(n) - V(k))*y(L) + b(L);

Ik = (V(k) - V(n))*y(L) + b(L)*V(k);

Snk = V(n)*conj(In)*basemva;

Skn = V(k)*conj(Ik)*basemva;

SL = Snk + Skn;

SLT = SLT + SL;

elseif tb(L)==n k = fb(L);

In = (V(n) - V(k))*y(L) + b(L)*V(n);

Ik = (V(k) - V(n))*y(L) + b(L);

Snk = V(n)*conj(In)*basemva;

Skn = V(k)*conj(Ik)*basemva;

SL = Snk + Skn;

SLT = SLT + SL;

else, end

if fb(L)==n || tb(L)==n

fprintf('%12g', k),

fprintf('%9.3f', real(Snk)), fprintf('%9.3f', imag(Snk))

fprintf('%9.3f', abs(Snk)),

fprintf('%9.3f', real(SL)),

if fb(L) ==n

fprintf('%9.3f \n', imag(SL))

else, fprintf('%9.3f\n', imag(SL))

end

else, end

end

end

SLT = SLT/2;

fprintf(' \n'), fprintf(' Total loss ')

fprintf('%9.3f', real(SLT)), fprintf('%9.3f\n', imag(SLT))

clear Ik In SL SLT Skn Snk

#exp-5 economic dispatch

clc

ng=input('enter the value of ng: ');

q=input('enter the value of load demand: ');

a=[300 200 201];

b=[5.3 5.5 5.6];

c=[0.004 0.006 0.009];

d=[450 350 225];

e=[200 150 100];

sum=0;

for i=1:1:ng

sum=sum+(b(i)/(2*c(i)));

end

sum=sum+q;

k=0;

for i=1:1:ng

k=k+(1/(2*(c(i))));

end

lambda =sum/k;

sn=0;

kn=0;

for i=1:1:ng

count=0;

pg(i)=(lambda-b(i))/(2*c(i));

if pg(i)>d(i)

rs(i)=d(i);

qnew=q-rs(i);

fprintf('the power at %d is %f\n',i,rs(i));

count=count+1;

end

end

for i=1:1:ng

if (pg(i)>d(i) || pg(i)<e(i))

for j=1:1:ng

if j~=i

sn=sn+(b(j)/(2*c(j)));

kn=kn+(1/(2*c(j)));

end

end

sn=sn+qnew;

lambda1=sn/kn;

disp(lambda1);

for (k=1:1:ng)

if (k~=i)

gs(k)=(lambda1-b(k))/(2*c(k));

fprintf('the final power at %d is %f\n',k,gs(k));

end

if (k==i)

fprintf('the final power at %d is %f\n',i,rs(i));

end

end

end

end

t=0;

for (i=1:1:ng)

h(i)=(a(i)+(b(i)*pg(i))+(c(i)*pg(i)*pg(i)));

fprintf('the total cost at %d is %f Rs/hr\n',i,h(i));

t=t+h(i);

end

fprintf('the fuel cost at id is %f Rs/hr\n',t);

#exp-6 fault analysis

clc

disp('creation of positive sequence matrix');

impedance;

Z1=Znew;

%disp('creation of negative sequence matrix');

Z2=Z1;

disp('creation of zero sequence matrix');

impedance;

Z0=Znew;

Z1

Z2

Z0

f=input('enter type of fault 1-symmetrical 2-LG fault');

if (f==1)

Zf=input('enter the fault impedance');

n=input('enter bus no at which fault occurs');

Zeq=Z1(n,n)+complex(0,Zf);

If=1/(abs(Zeq));

fpritnf('the fault current is %d',If);

else if (f==2)

Zf=input('enter the fault impedance');

n=input('enter bus no at which fault occurs');

Zeq=Z1(n,n)+Z2(n,n)+Z0(n,n)+3*complex(0,Zf);

If=3/(abs(Zeq));

fpritnf('the fault current is %d',If);

end

end

Impedance;

n=input('enter no.of nodes including reference node: ');

link=input('enter bo of co-trees or links');

z=(zeros(n-1,n-1));

del=(zeros(n-1,1));

for i=1:1:n-1

a=input('enter to which node you want to attach the branch: ');

for j=1:1:n-1

if a==0

c=input('enter resistance');

d=input('enter reactance');

z(i,i)=complex(c,d);

break;

else

if i==j

c=input('enter resistance');

d=input('enter reactance');

z(i,j)=z(a,a)+complex(c,d);

else

z(i,j)=z(a,j);

z(j,i)=z(j,a);

end

end

end

end

Zold=z;

for k=1:1:link

disp('enter x,y node where you want to add co-tree: ');

x=input("x=");

y=input("y=");

c=input('enter resistance');

d=input('enter reactance');

z(n,n)=complex(c,d)+z(x,x)+z(y,y)-2*z(x,y);

for j=1:1:n-1

z(n,j)=z(y,j)-z(x,j);

z(j,n)=z(j,y)-z(j,x);

del(j,1)=z(j,n);

end

m=(del*del.')/z(n,n);

Znew=Zold-m;

Zold=Znew;

z=Znew;

end

%disp('the Zbus matrix is');

Znew;