1. To find ABCD constants of equivalent PI transmission line and hence
efficiency and regulation ALGORITHM: 1. Read length, impedance/km(zkm), admittance/km of transmission line(ykm) 2. Enter choice 1-for short line 2-for medium line 3-for long line 3. If choice=1 Read length Set B = Z A = D = 1 C = 0 If choice=2 Read length A=D=1 + yz/2 B = Z C = y (1 + yz/4) If choice=3 Read length Zeq zc = z/y = z*y*length A = D = Cosh () Yeq Yeq B = Zc Sinh () C = 1/Zc Sinh () Find generalized equivalent circuit elements of network zeq = z*length*sinh ()/ yeq = y*length/2*tanh (/2)/ (/2) Display Zeq & Yeq 4. Display A, B, C, & D constants 5. Read choice 1 – To read V r, I r and compute V s , I s 2 – To read V s, I s and compute V r, I r 6. If choice = 1 Find V s = A V r + B I r I s = C V r + D I r where V r & I r are in per ph ase Display V s & I s. if choice = 2 Find V r = D Vs – B Is I r = -C Vs+ A Is where Vs & Is are in per phase Display V r & I r. 7. Find (i) VS SV & Is SI (ii) Sending end Power factor Pf-send = cos (SV -SI) (iii) Sending end Power, PS = 3 VS x I s x Pf-send 8. Find (i) V r rv & I r rI (ii) Receiving end power factor Pf-rec = cos (rV -rI) (iii) Receiving end power, Pr=3 V r x I r x Pf-rec 9. Find efficiency, = Pr x 100 Ps 10.
Find voltage regulation VR = Vs
-
V r
X 100 V r 11. 12.
Display (i) efficiency, (ii) Voltage regulation, VR Stop
PROGRAM: % ABCD of eqt. PI Network zkm = 0.2+0.408i ; ykm= 0+3.14e-6i; k1=input ('\n Enter 1-for short line 2-for medium line 3-for long line switch k1 case 1, length = input ('\n Enter length of short line '); Z=zkm*length; Y=ykm*length; A=1; B=Z; C=0; D=1; case 2, length = input ('\n Enter length of medium line '); Z=zkm*length; Y=ykm*length; A=1+Y*Z/2; B=Z; C=Y*(1+Y*Z/4); D=A; case 3, length = input ('\n Enter length of long line '); zc=sqrt(zkm/ykm); gam=sqrt(zkm*ykm)*length; A=cosh (gam); D=A; B=zc*sinh(gam); C=1/zc*sinh(gam); fprintf ('\n The equivalent PI circuit constants :'); zeq = zkm*length*sinh (gam)/gam; yeq = ykm*length/2*tanh (gam/2) / (gam/2); fprintf ('\n Zeq = %15.4f %+15.4fi',real(zeq),imag(zeq)); fprintf ('\n Yeq/2= %15.4f %+15.4fi',real(yeq),imag(yeq)); otherwise disp('wrong choice of tr.line'); end fprintf('\nA,B,C and D constants : \n'); fprintf('----------------------------------------'); fprintf('\nA = %15.4f %+15.4fi',real(A),imag(A)); fprintf('\nB = %15.4f %+15.4fi',real(B),imag(B)); fprintf('\nC = %15.4f %+15.4fi',real(C),imag(C)); fprintf('\nD = %15.4f %+15.4fi',real(D),imag(D)); fprintf('\n The product AD-BC=%f',A*D-B*C); k2=input('\n Enter 1 - To read Vr, Ir and compute Vs , Is \n 2 To read Vs, Is and compute Vr, Ir '); switch k2, case 1, vr=input('enter Vr/line '); ir=input('enter Ir/phase '); vr=vr*1e3/sqrt(3); vs = (A*vr+B*ir)/1e3; is = C*vr+D*ir; fprintf('\nSending end Voltage/ph=%f %+fi KV',real(vs),imag(vs)); fprintf('\nSending end Current/ph=%f %+fi AMP',real(is),imag(is)); vs=vs*1e3; case 2, vs=input('enter Vs/line '); is=input('enter Is/phase '); vs=vs*1e3/sqrt(3.0); vr=(D*vs-B*is)/1e3; ir= -C*vs+A*is; fprintf('\nReceiving end Voltage/ph=%f %+fi KV',real(vr),imag(vr));
');
fprintf('\nReceiving end Current/ph=%f %+fi AMP',real(ir),imag(ir)); vr=vr*1e3; otherwise disp('wrong choice'); end rec_pow=3*real(vr*conj(ir))/1e6; send_pow=3*real(vs*conj(is))/1e6; eff=rec_pow/send_pow*100; reg=(abs(vs)/abs(A)-abs(vr))/abs(vr)*100; fprintf('\n Receiving end power=%.2f KVA',rec_pow); fprintf('\n Sending end power=%.2f KVA',send_pow); fprintf('\n Efficiency=%.2f %%',eff); fprintf('\n Voltage Regulation=%.2f%%',reg); 2. To find ABCD constants of equivalent T transmission line and hence efficiency and regulation ALGORITHM: 1. Read length, impedance/km, admittance/km of transmission line 2. Enter choice 1-for short line 2-for medium line 3-for long line 3. If choice=1 Read length Set B = Z Zeq A = D = 1 C = 0 If choice=2 Read length Set A=D=1 + YZ/2 B = Z* (1 + YZ/4) C =YZ If choice=3 Read length zc = z/y = (z*y)*length A = D = Cosh () B = zc Sinh () C = 1/zc Sinh () zeq=z*length/2*tanh(/2)/( /2); yeq=y*length*sinh()/; Display Zeq & Yeq 4. Display A, B, 5. Read choice 1 2 – To read V s, I s 6. If choice = 1 Find V
Zeq Yeq
C, & D constants – To read V r, I r and compute V s , I s and compute V r, I r s = A V r + B I r I s = C V r + D I r
where V r & I r are in per p
hase Display V s & I s. if choice = 2
7. 8.
Find V r = D Vs – B Is I r = -C Vs+ A Is where Vs & Is are in per phase Display V r & I r. Find (i) VS SV & Is SI (ii) Sending end Power factor Pf-send = Cos (SV -SI) (iii) Sending end Power, PS = 3 VS x I s x Pf-send Find (i) V r rv & I r rI (ii) Receiving end power factor Pf-rec = Cos (rV -rI)
9. 10.
(iii) Receiving end power, Pr=3 V r x I r Find efficiency, = Pr x 100 Ps Find voltage regulation VR = Vs - V r
x Pf-rec
X 100 V r 11. 12.
Display (i) efficiency, (ii) Voltage regulation, VR STOP
PROGRAM: %ABCD of eqt. T Network zkm =0.2+0.408i; ykm=0+3.14e-6i; k1=input('\n Enter 1-for short line 2-for medium line 3-for long line switch k1 case 1, length=input('\n Enter length of short line '); Z=zkm*length; Y=ykm*length; A=1;B=Z;C=0;D=1; case 2, length=input('\n Enter length of medium line '); Z=zkm*length; Y=ykm*length; A=1+Y*Z/2; B=Z; C=Y*(1+Y*Z/4); D=A; case 3, length=input('\n Enter length of long line '); zc=sqrt(zkm/ykm); gam=sqrt(zkm*ykm)*length; A=cosh(gam); D=A; B=zc*sinh(gam); C=1/zc*sinh(gam); otherwise disp('wrong choice of tr.line'); end fprintf('\nA,B,C and D constants : \n'); fprintf('----------------------------------------------------'); fprintf('\nA = %15.4f %+15.4fi',real(A),imag(A)); fprintf('\nB = %15.4f %+15.4fi',real(B),imag(B)); fprintf('\nC = %15.4f %+15.4fi',real(C),imag(C)); fprintf('\nD = %15.4f %+15.4fi',real(D),imag(D)); k2=input('\n Enter 1 - To read Vr, Ir and compute Vs , Is \n 2 To read Vs, Is and compute Vr, Ir '); switch k2, case 1, vr=input('enter Vr/phase '); ir=input('enter Ir/phase '); vr=vr*1e3/sqrt(3); vs=(A*vr+B*ir)/1e3; is=C*vr+D*ir; fprintf('\nSending end Voltage/ph=%f %+fi KV',real(vs),imag(vs)); fprintf('\nSending end Current/ph=%f %+fi AMP',real(is),imag(is)); vs=vs*1e3; case 2, vs=input('enter Vs/phase '); is=input('enter Is/phase '); vs=vs*1e3/sqrt(3.0);
');
vr=(D*vs-B*is)/1e3; ir=-C*vs+D*is; fprintf('\nReceiving end Voltage/ph=%f %+fi KV',real(vr),imag(vr)); fprintf('\nReceiving end Current/ph=%f %+fi AMP',real(ir),imag(ir)); vr=vr*1e3; end rec_pow=3*real(vr*conj(ir))/1e6; send_pow=3*real(vs*conj(is))/1e6; eff=rec_pow/send_pow*100; reg=(abs(vs)/abs(A)-abs(vr))/abs(vr)*100; fprintf('\n Receiving end power=%.2f KVA',rec_pow); fprintf('\n Sending end power=%.2f KVA',send_pow); fprintf('\n Efficiency=%.2f %%',eff); fprintf('\n Voltage Regulation=%.2f%%\n',reg); 3. Power Angle Characteristics of Cylindrical Rotor Generator ALGORITHM: 1. Read x, V and Ia 2. E = V + j x Ia Pm = E V / Xd 3. Plot power angle characteristics from delta 0 to pi Indicate Pmax and corresponding delta E - V 4. Regulation R = ---------x 100 V Display R 4. Stop PROGRAM: %power angle characteristics cylindrical rotor syn.m/c % enter the values of x,v,Ia; x=0.3; v=1.0; ia=0.5000-0.375i; e=v+x*ia*i; pm=abs(e)*abs(v)/x; del=0:0.1:pi; plot(del,pm*sin(del),'r'); title('power angle characteristics cylindrical rotor syn.m/c'); xlabel('Delta'); ylabel('power ,p'); legend('pm*sin(delta)'); del=0:0.01:pi/2; line(del,pm); y=0:0.01:pm; line(pi/2,y); reg=(abs(e)-abs(v))/abs(v)*100; fprintf('\n the excitation voltage=%f%%\n',e); fprintf('\n the peak power=%f%%\n',pm); fprintf('\n the voltage regulation=%f%%\n',reg); 4. Power Angle characteristics of salient rotor generator ALGORITHM: 1. Read Xd, Xq, Xd1, V, P & Pf 2. Ө = Cos-1(pf) I = P / (V x pf) Ia = I cosӨ - j I SinӨ delta = tan((XE) Ia CosӨ ) / (IV) + Xq Ia SinӨ)) E = V Cos(delta) + Xd / Ia Sin (delta + Ө) P = E V / Xd
Pr = V 2 (Xd – Xq) / (2 x d x q) Pe = P Sm(delta) + Pr Sin ( 2 delta) 3. Plot P, Pr, Pe from delta = 0 to Pi Find max Pe and the corresponding delta and indicate the same E - V 4. Regulation R = --------------- x 100 V Display E , maxPe, delta ,Reluctance power (Pr) & R 5. Stop PROGRAM: %power angle characteristics salient pole rotor syn.m/c % enter the values of xd,xq,xdd,v,Ia; xd=1.0 ; xq=0.6; v=1.0; ia=0.5000-0.375i; theta=abs(angle(ia)); del=atan((xq*abs(ia)*cos(theta))/(abs(v)+xq*abs(ia)*sin(theta))); e=abs(v)*cos(del)+xd*abs(ia)*sin(del+theta); p=e*abs(v)/xd; pr=abs(v)^2*(xd-xq)/(2*xd*xq); del=0:0.01:pi; line(del,0);hold on; pe=p*sin(del)+pr*sin(2*del); plot(del,p*sin(del), 'r',del,pr*sin(2*del),'m',del,p*sin(del)+pr*sin(2*del),'b'); title('power angle characteristics salient pole syn.m/c'); xlabel('Delta'); ylabel('power ,p'); legend('p*sin(delta)','pr*sin(2delta)','p*sin(delta)+pr*sin(2delta)'); [pmax,k]=max(pe); fprintf('\n maximum power=%f and the corresponding angle=%f',pmax,del(k)*180/pi) del=0:0.01:del(k); line(del,pmax); y=0:0.01:pmax; line(del(k),y); reg=(abs(e)-abs(v))/abs(v)*100; text(-0.4,pmax,'pmax'); text(del(k),-0.1,'deltamax'); fprintf('\n the excitation voltage =%f pu',abs(e)); fprintf('\n the voltage regulation=%f%%',reg); fprintf('\n Reluctance power= %f p.u\n',pr); 5. Swing curve under sustained fault ALGORITHM: 1. Read Input power (P s), voltages (e & v), M, External reactance (xe), line reactance (xl) 2. Pre fault reactance,x1=xe+xL/2 3. Read choice 1-for fault at the beginning 2-for fault at the middle if choice=1 x2=infinity if choice=2 x2= (xe*xL+xe*xL/2+xl*xL/2)/(xL/2); 4. Set dt=0.001 r r=180/ (required to convert deg into radius etc) f=(dt)2 / M Iteration Count, it=1 Time, t(1)=0 =0 5. Find Pm1 = e1*e2/X1 0 = sin-1(Ps/Pm1) Pm2 = Pm=e1*e2/X2 Pe = Pm Sin(0)
Pa = (Ps-Pe)/2 Repeat step no. i to step no. vi from time=0sec to 1 sec in steps of 0.0
6. 5 sec I. Find ft = f*Pa/rr II. = +ft III. Display time (t), Pmax (Pm), in degrees. IV. Increment iteration count it=it + 1 V. Find t(it)=t(it-1) + dt VI. (it) = (it-1)+ VII. Set Pm=Pm2 VIII. Find Pe=Pm Sin IX. Pa=Ps-Pe 7. Plot swing curve 8. Find critical clearing angle c=cos-1(ps(m-o)+pm3*cos(m)-pm2*cos(o)) / (pm3-pm2) Where pm3=e*v/x3, x3=xe+x1, m=pi-o. 9. Find the critical clearing time(Tc), corresponding to c and dislay Tc 10. Stop PROGRAM: %swing curve for sustained fault and critical clearing angle & time. %ps=mech.power input xe=xg+xt, x1=reactance before fault x2=reactance after %fault clc; clear all; % ps = 0.9; e1 = 1.1; v = 1.0; m = 0.00028; x1 = 0.45; x2 = 1.25; x3 = 0.55; data = input( ‘ps = , e1= , v= , m= , x1= , x2= , x3= in matrix form’); ps = data(1); e1 = data(2); v = data(3); m = data(4); x1 = data(5); x2 = data(6); x3= data(7); rr = 180/pi; dt = 0.05; rr=180/pi; f=dt^2/m; it=1; t(it)=0; deld=0; pm1=e1*v/x1; del(it)=asin(ps/pm1); pm2=e1*v/x2; pm=pm2; pe=pm*sin(del(it)); pa=(ps-pe)/2; fprintf( \n SUSTAINED FAULT ); fprintf( \n---------------------------- ); fprintf( \n TIME PMAX DELTA ); fprintf( \n---------------------------- ); t1=0; while(t(it)<=1.0) ft=f*pa/rr; deld=deld+ft; if(t1-t(it)<0.0001) fprintf( \n %5.3f %5.2f %5.2f ,t(it),pm,del(it)* rr); t1=t1+0.05; end it=it+1; t(it)=t(it-1)+dt; del(it)=del(it-1)+deld; pm=pm2; pe=pm*sin(del(it));
pa=ps-pe; end plot(t,del, r ); title( swing curve ); xlabel( Delta ); ylabel( power ,p ); pm3=e1*v/x3; delm=pi-del(1); cdc=(ps*(delm-del(1))+pm3*cos(delm)-pm2*cos(del(1)))/(pm3-pm2); delc=acos(cdc); fprintf( \n critical clearing angle-%f ,delc*rr); it=1; while(t(it)<1.0) if (del(it)>=delc) break; end it=it+1; end fprintf( \n critical clearing time-%f ,t(it)); 6. Swing curve with fault being cleared ALGORITHM: 1. Read Input power (P s), voltages (e & v), M, External reactance(xe),line reactance (xl) pre fault reactance,x1=xe+xL/2 2. Read choice 1-for fault at the beginning 2-for fault at the middle if choice=1 x2=infinity if choice=2 x2= (xe*xL+xe*xL/2+xL*xL/2)/(xL/2); 3. Set x3=xe+xL dt=0.05 4. Read clearing time ct in sec’s K=ct/dt (integer division) r=ct-k*dt If r=0 then Print ‘Fault is cleared at the beginning of an interval’ else Print ‘Fault is cleared in the middle of an interval’ 5. Set rr=180/ f=dt2/M Iteration Count, it=0 Time, t=0 6. Find Pm1=e1*e2/x1 0=Sin-1 Ps/Pm1 Pm2=e1*e2/x2 Pm3=e1*e2/x3 Pm=Pm2 Pe=Pm * Sin (0) Pa=Ps-Pe 7. Repeat step no. i to step no.vii from time t=0 to 1 sec in steps of 0.05 sec i. Find ft=f * Pa/rr = +ft ii. Display time (t), Pmax (Pm), in degrees. iii. Increment iteration count it=it+1 iv. Find t(it)=t(it-1) + dt i. (it)= (it-1)+ v. If r=0 then
i. If(time < clearing time) then 1. Pm=Pm2 ii. else if(time=clearing time) then 1. Pm=(Pm2+Pm3)/2 b. If(r! =0) then i. If (time < clearing time) then 1. Pm=Pm2 ii. else 1. Pm=Pm3 vi. Find Pe=Pm Sin vii. Pa=Ps-Pe 8. Plot swing curve 9. If decreases at any point then the system is stable otherwise unstable. Find whether the system is stable or not and display the same. 10. Stop PROGRAM: clc; clear all; %program for swing curve when the fault is cleared % ps=0.9; e=1.1; v=1.0; m=0.00028; % x1=0.45; x2=1.25; x3=0.55; data = input( ‘ps = , e1= , v= , m= , x1= , x2= , x3= in matrix form’); ps = data(1); e1 = data(2); v = data(3); m = data(4); x1 = data(5); x2 = data(6); x3= data(7); dt=0.05; ct=input( \n enter clearing time in secs ); k=ct/dt; r=ct-floor(k)*dt; if(r==0) fprintf( FAULT IS CLEARED AT THE BEGINNING OF AN INTERVEL ); else fprintf( FAULT IS CLEARED AT THE MIDDLE OF AN INTERVEL ); end ct = ct/50; rr=180/pi; f=dt*dt/m; deld=0; it=1; t(it)=0; pm1=e*v/x1; del(it)=asin(ps/pm1); pm2=e*v/x2; pm3=e*v/x3; pm=pm2; pe=pm*sin(del(it)); pa=(ps-pe)/2; fprintf( \n TIME PMAX DELTA ); while(t(it)<=1.05) ft=f*pa/rr; deld=deld+ft; fprintf( \n %5.2f %5.2f %5.2f ,t(it),pm,del(it)* rr); it=it+1; t(it)=t(it-1)+dt; del(it)=del(it-1)+deld; if(r==0) if(t(it)
else pm=pm3; end end if(r~=0) if(t(it)
R = z(:,3); X = z(:,4); nbr = length(fbus); nbus = max(max(fbus),max(tbus)); Z = R+j*X; y = ones(nbr,1)./Z; Ybus = zeros(nbus,nbus); for k = 1: nbr if fbus(k)>0 & tbus(k)>0 Ybus(fbus(k),tbus(k)) = -y(k); Ybus(tbus(k),fbus(k)) = Ybus(fbus(k),tbus(k)); else, end end Ybus; for n = 1:nbus for k = 1:nbr if fbus(k)==n tbus(k)==n; Ybus(n,n) = Ybus(n,n)+y(k); else, end end end Ybus 8. Formation of Ybus using singular transformation method with half line charging (or shunt admittance)
ALGORITHM: 1. Read no. of buses( nbus) and no. of lines( nline). 2. Initialize shunt admittance matrix to zeros. 3. Repeat step no i. to step no iv for I=1 to nline. i. Read line data i.e., frombusno, tobus no, Zseries, half line charging. ii. Shunt admittance [from bus]=shunt admittance [from bus]+half line chargi ng. iii. Shunt admittance [to bus]=shunt admittance [to bus]+half line charging. iv. Series admittance [I]=1/series impedance [I]. 4. Form the incidence matrix A (i) Initialize all the elements of matrix A of size [(nbus + nline) x nbus] to zeros (ii) Form A = Where AI is identity matrix of size (nbus x nbus) AA is a matrix of size (nline x nbus) and elements of AA = [aij] Where aij = 0 if ith branch is not incident on jth bus. aij = 1 if ith branch is incident to and oriented away from jth bus. aij = -1 if ith branch is incident to and oriented away from jth bus. 5. Form Yprimitive matrix Initialize all elements of matrix Ypr of size (nbus+nline x nbus+nline) to zeros Form Ypr
=
Where Y is a matrix of size (nbus x nbus) with diagonal elements Yii = total shunt admittance at bus i and all other elements being zeros.YY is a matrix of size (nbus x nline) with all elements being zeros.Y1 is matrix of size (nline x nline) with diagonal elements of Y1ii =Yseries of line i=1/Zseries of line i. YY1 is matrix of size (nline x nbus) with all elements being zeros.
6. 7. 8. 9. 10. 11.
Display A matrix & Ypr matrix. Form Transpose of matrix A i.e., AT. Form ATY = AT x Ypr. Form Ybus = ATY x A and display. Form Zbus = [Ybus]-1 and display. Stop
PROGRAM: %FORMATION OF ybus USING SINGULAR TRANSFORMATION METHOD %WITHOUT MUTUAL COUPLING: % p q Z hlcY(ADM) z=[
1 2 0.02+0.06i 0.03i 1 3 0.08+0.24i 0.025i 2 3 0.06+0.18i 0.02i 2 4 0.06+0.18i 0.02i 2 5 0.04+0.12i 0.015i 3 4 0.01+0.03i 0.01i 4 5 0.08+0.24i 0.025i fb=z(:,1); tb=z(:,2); Z=z(:,3); hlcy=z(:,4); y=1./Z; nbus=max(max(fb),max(tb)); Y=zeros(nbus); nline=length(fb); nlb=nline+nbus; A=zeros(nlb,nbus); for k=1:nbus A(k,k)=1; end for k=1:nline A(nbus+k,fb(k))=1; A(nbus+k,tb(k))=-1; end sh=zeros(nbus); for k=1:nline sh(fb(k))=sh(fb(k))+hlcy(k); sh(tb(k))=sh(tb(k))+hlcy(k); end ypr=zeros(nlb,nlb); for k=1:nbus ypr(k,k)=sh(k); end for k=1:nline ypr(nbus+k,nbus+k)=y(k); end format short; Ybus=A'*ypr*A
9.
];
Formation of Ybus using singular Transformation method with mutual coupling.
ALGORITHM: 1. Read no. of buses( nbus) and no. of lines (nline). 2. Road line data i.e., Frombus no, tobus no, Zseries, mutual line no, Zmutual at all lines. 3. Form bus incidence matrix A (i) Initialize all the elements of A of size (nline x nline) to zeros (ii) Form A = [aij] Where aij = 0 if ith branch is not incident on jth bus. aij = 1 if ith branch is incident to and oriented away from jth bus. aij = -1 if ith branch is incident to and oriented away from jth bus. 4. Form Zprimitive matrix, Zpr (iii) Initialize all the elements of Zpr of size (nline x nline) to zeros (iv) The diagonal elements of Zpr are assigned with respective Zseries i.e., Zpr (i, i) = Zseries of ith line (v) The off diagonal element is zero if the line i has no mutual coupling with line j (vi) The off diagonal element Zpr (i, j) = Zmutual if the line i has mutual coupling with line j 5. Form Y primitive matrix Ypr = [Zpr]-1 6. Display A matrix & Ypr matrix. 7. Form Transpose of matrix A i.e., AT. 8. Form ATY = AT x Ypr. 9. Form Ybus = ATY x A and display. 10. Form Zbus = [Ybus]-1 and display. 11. Stop PROGRAM: %FORMATION OF ybus USING SINGULAR TRANSFORMATON METHOD WITH MUTUAL COUPLING & WITHOUT LINE CHARGING: % p q Z mno mutual(imp) z=[ 0 1 0.6i 0 0 0 2 0.5i 1 0.1i 2 3 0.5i 0 0 0 1 0.4i 1 0.2i 1 3 0.2i 0 0 ]; p=z(:,1); q=z(:,2); Z=z(:,3); mno=z(:,4); zmc=z(:,5); nbus=max(max(p),max(q)); Y=zeros(nbus); nline=length(p); A=zeros(nline,nbus); for k=1:nline, if(q(k) = = 0) A(k,p(k))=1; elseif(p(k) = = 0) A(k,q(k)) = -1; end if(p(k)~=0 & q(k)~=0) A(k,p(k))=1; A(k,q(k))=-1; end end zpr=zeros(nline,nline); for k=1:nline zpr(k,k)=Z(k); if(mno(k))
zpr(k,mno(k))=zmc(k); zpr(mno(k),k)=zmc(k); end end ypr=inv(zpr); format short; Ybus=A'*ypr*A 10.
Bus cuurents, Bus power and Line Flows
ALGORITHM: 1. Read the following line data at each line from bus no(p), to bus no(q), line impedance(Zpq) and half line charging admittance(Ychpq) Read voltages at all buses 2. Convert line impedance to admittance of each line Ypq 3. Repeat step i to ix for K=1 to no. of lines i. P = from bus no q = tobusno ii, Line current, Ipq = (Vp-Vq)Ypq + Vq x Ychpq iii. Line loss = SLpq =VpIpq* iv. Generation at bus p = SGp = SGp + S2P2 v. Iqp = (Vq-Vp)Ypq + Vq x Ychqp vi. Iq = Iq + Iqp vii. SLqp + Vq Iqp* viii. Loss (k) = Loss + SLpq + SLqp ix. Total loss = Total loss + loss(k) 5. Display current at each bus Display line flows and line loss at each line Display generations at each bus Display the total losses 6. Stop
PROGRAM: % p q
z
hlc(ADM)
z= [
5 4 0.02+0.06i 0.03i 5 1 0.08+0.24i 0.025i 4 1 0.06+0.18i 0.02i 4 2 0.06+0.18i 0.02i 4 3 0.04+0.12i 0.015i 1 2 0.01+0.03i 0.01i 2 3 0.08+0.24i 0.025i ]; v(1)=1+0.2i; v(2)=1.01-0.5i; v(3)=1.05-0.3i; fb=z(:,1); tb=z(:,2); Z=z(:,3); ysh=z(:,4); y=1./Z; nbus=max(max(fb),max(tb)); nline=length(fb); II=zeros(nbus,1);SG=zeros(nbus,1); totloss=0; for k=1:nline p=fb(k);q=tb(k);
v(4)=1.01-0.5i;
v(5)=1.0-0i;
I(p,q)=(v(p)-v(q))*y(k)+v(p)*ysh(k); II(p)=II(p)+I(p,q); SL(p,q)=v(p)*conj(I(p,q)); SG(p)=SG(p)+SL(p,q); I(q,p)=(v(q)-v(p))*y(k)+v(q)*ysh(k); II(q)=II(q)+I(q,p); SL(q,p)=v(q)*conj(I(q,p)); SG(q)=SG(q)+SL(q,p); loss(k)=SL(p,q)+SL(q,p); totloss=totloss+loss(k); end f=fopen('lf.out','w'); fprintf(f,'bus currents'); fprintf(f,'\nBus no Current'); for k=1:nbus fprintf(f,'\n %d %10.4f %+10.4fi',k,real(II(k)),imag(II(k))); end fprintf(f,'\nLine flows'); fprintf(f,'\nFrom bus To bus Lineflows Line Loss'); for k=1:nline p=fb(k);q=tb(k); fprintf(f,'\n %d %d %10.4f%+10.4fi %10.4f%+10.4fi',p,q,real(SL(p,q)),imag(SL(p,q)),real(loss(k)),imag(loss(k))); fprintf(f,'\n %d %d %10.4f%+10.4fi %10.4f%+10.4fi',q,p,real(SL(q,p)),imag(SL(q,p)),real(loss(k)),imag(loss(k))); end fprintf(f,'\nBus Generations'); fprintf(f,'\nBus no Generation'); for k=1:nbus fprintf(f,'\n%d %10.4f%+10.4fi',k,real(SG(k)),imag(SG(k))); end fprintf(f,'\n\nTotal Losses=%10.4f%+10.4fi',real(totloss),imag(totloss)); fclose(f);
11. Load flow analysis using Gauss-Seidal method. ALGORITHM: 1. Read Ybus, no.of buses (nbus) (Assume bus1 as slack bus) No.of pq buses (Npq), No.of pv buses (Npv) Voltage at slack bus Real power (P) and reactive power (Q) at all pq buses Real power (P) and specified voltage (Vsp) and Qlimits (Qmin and Qmax) at PV buses. Convergance critirian ( ). 2. Set iteration count k=1. 3. Calculate voltages at all the pq buses( p=2 to Npq+1), using 4. To calculate voltages at all the pv buses( p=Npq+2 to nbus) Find QCp= -Im( ) whare If QCp is within Qmin and Qmax, treat the bus p as pv bus and calculate where at an angle p corresponding to previous iteration Else if QCp violates Qlimits, treat bus p as pq bus and if QCP< Qmin , set Qp=Qmin else set Qp=Qmax and find where corresponding to previous iteration
5. 6. 7. 8.
Display voltages at all buses and QCp at pv buses Determine the largest absolute value of . If goto next step otherwise set k=k+1and goto step no.3. Stop
PROGRAM: % enter Y-bus of the system Y=[ 20-50i -10+20i -10+30i -10+20i 26-52i -16+32i -10+30i -16+32i 26-62i ]; %enter no.of buses,slack bus no,no.of pq buses,no of pv buses nbus=3; sbno=1; npq=1; npv=1; % assume bus no 1 as slack bus ,next npq buses as pq buses and the remaining as pv buses %enter slack bus voltage v(1)=1.05+0i; %assume voltages at all other buses as 1+j0 v(2)=1; v(3)=1; %enter p & q at all pq buses P(2)= -4.0; Q(2)= -2.5; %enter p ,vsp and q limits at pv buses P(3)=2.0; vsp(3)=1.04; Qmin(3)=0.1; Qmax(3)=1.0; %enter accuracy of convergence acc=0.001; for it=1:5 %to find voltages at pq buses for p=2:npq+1, v1(p)=v(p); pq=(P(p)-Q(p)*i)/conj(v(p)); ypq=0; for q=1:nbus if(p= =q) continue; end ypq=ypq+Y(p,q)*v(q); end v(p)=(pq-ypq)/Y(p,p); end %to find voltages at pv buses for p=npq+2:nbus, v1(p)=v(p); s=0; for q=1:nbus, if(p~=q) s=s+Y(p,q)*v(q); else vp=v(q); ang=angle(v(q)); v(q)=complex(vsp(q)*cos(ang),vsp(q)*sin(ang)); s=s+Y(p,q)*v(q); end end Qc(p)=-1*imag(conj(v(p))*s); if(Qc(p)>=Qmin(p)& Qc(p)<=Qmax(p)) pq=(P(p)-Qc(p)*i)/conj(v(p)); ypq=0;
for q=1:nbus if(p= =q) continue; end ypq=ypq+Y(p,q)*v(q); end v(p)=(pq-ypq)/Y(p,p); ang=angle(v(p)); v(p)=complex(vsp(p)*cos(ang),vsp(p)*sin(ang)); else if(Qc(p)
12. Formation of Jacobin Matrix in polar form without PV Bus ALGORITHM: 1. Read (i) Y bus (ii) No. of buses,n (iii) Voltages at all buses 2. . Repeat for p =1 to (n-1) Өpq = Өp - Өq where Өp = ∟Vp & Өq = ∟Vq Gpq + real (Ypq), Bpq = image (Ypa) Pcal(p) = Vp
(Gpq + Bpq SinӨpq) Vq)
Qcal(p) = Vp
(GpqSm Өpq - BpqCos Өpq)Vq
3.
Off diagonal of J1 = VpVq(Gpq SinӨpq-BpqCosӨpq) = off diagonal elements of J4 Off diagonal of J2 = VpVq(Gpq CosӨpq + Bpq Sin Өpq) = - Off diagonal elements of J3 Diagonal of J1 = -Qcal (p) – Bpp Vp2 J2 = Pcal (p) + Gpp Vp2 J3 = Pcal (p) – Gpp Vp 2 J4 = Qcal (p) – Bpp Vp2 4. Combine J1, J2, J3 & J4 into a single matrix J 5. Display J 6. Stop PROGRAM: % enter Y-bus of the system Y=[
20-50i -10+20i -10+30i -10+20i 26-52i -16+32i -10+30i -16+32i 26-62i ]; %enter no.of buses nbus=3; % assume last bus as slack bus ,and the remaining as pq buses %assume voltages at all other buses as 1+j0 v(1)=1+0i; v(2)=1+0i; %enter slack bus voltage v(3)=1.05+0i; %to find pcal and qcal for p=1:nbus-1 vp=abs(v(p)); pc(p)=0;qc(p)=0; for q=1:nbus vq=abs(v(q)); thpq=angle(v(p))-angle(v(q)); Gpq=real(Y(p,q));Bpq=imag(Y(p,q)); pc(p)=pc(p)+(Gpq*cos(thpq)+Bpq*sin(thpq))*vq; qc(p)=qc(p)+(Gpq*sin(thpq)-Bpq*cos(thpq))*vq; end pcal(p)=pc(p)*vp;qcal(p)=qc(p)*vp; end %formation of J1 J2 J3 & J4 for p=1:nbus-1 vp=abs(v(p));Gpp=real(Y(p,p));Bpp=imag(Y(p,p)); for q=1:nbus-1 vq=abs(v(q)); thpq=angle(v(p))-angle(v(q)); Gpq=real(Y(p,q));Bpq=imag(Y(p,q)); if(p~=q) J1(p,q)=vp*vq*(Gpq*sin(thpq)-Bpq*cos(thpq)); J2(p,q)=vp*vq*(Gpq*cos(thpq)+Bpq*sin(thpq)); J3(p,q)=-J2(p,q); J4(p,q)=J1(p,q); end end J1(p,p)=-qcal(p)-Bpp*vp^2; J2(p,p)=pcal(p)+Gpp*vp^2; J3(p,p)=pcal(p)-Gpp*vp^2; J4(p,p)=qcal(p)-Bpp*vp^2; end
%to combine J1,J2,J3&J4 into single matrix J J = [J1 J2; J3 J4]; fprintf('The Jacobian matrix is \n'); disp(J);
13. Z bus by Z- BUS Building Algorithm ALGORITHM: Type 1 : Line added between New bus to Reference bus Type 2 : Line added between New bus to Old bus Type 3 : Line added between Old bus to reference bus Type 4 : Line added between two Old buses Note : [ Treat BUS – 1 as Reference bus Form the input data from single line diagram. Start with type 1 only. Select starting bus (SB) as either Reference bus or existing bus ( for type 1 or 2 respectively ) and Select ending bus (EB) as new bus .] Initialize zbus with the size 1 X 1 size =nsize=rows=column=1; Read line data , depending on the type of the line specified, form the Z-bus as follows Type - 1 : Increase the size of zbus by 1 zbus( size,size) = zline; % diagonal element all the other off diagonal elements (added row & column) being zero Type – 2: Increase the size of zbus by 1 zbus( size,size) = zbus ( sb , sb ) + zline; -----diagonal element for off diagonal elements Copy contents of the row corresponding to sb to new row and to new column Type - 3 : size remains the same zsum = zbus( sb , sb ) + zline(i); repeat the steps i to ii for k= 2 to nsize step i : repeat the step ii for j = 2 to nsize step ii : zmod(k,j) = -1 * zbus(k,m) * zbus( m, j) / zsum; then the modified z bus is zbus=zbus+zmod; Type - 4 : size remains the same Calculate zsum = + zline (i) + zbus ( sb , sb )+zbus ( eb , eb )-2 * zbus( sb , eb ) repeat the steps i to ii for k= 2 to nsize step i : repeat the step ii for j = 2 to nsize step ii : zmod (k , j )=-1* ( zbus( k , sb ) - zbus(k , eb ) ) * ( zbus (sb , j)-zbus( eb , j )) / zsum the modified z bus is zbus = zbus + zmod PROGRAM: % Z-BUS --- using Z-Bus Building Algorithm % data for line = line no , type, SB, EB, line_impedance zdata=[ 1 1 1 2 1.0 2 2 2 3 1.0
3 2 2 3 1.0 4 3 3 1 1.0 5 3 4 1 1.0]; line=zdata(:,1); type=zdata(:,2); sb=zdata(:,3); eb=zdata(:,4); zline=zdata(:,5); nb=max(max(sb),max(eb)); nl= length(line); %zbus=zeros(nb); % formation of Z-BUS nsize=1; for i= 1:nl m=sb(i); p=eb(i); % type 1 if (type(i)==1) nsize=nsize+1; zbus(nsize,nsize)=zline(i); end % type 2 if (type(i)==2) nsize=nsize+1; for j=2:(nsize-1) zbus(nsize,j)=zbus(m,j); zbus(j,nsize)=zbus(j,m); end zbus(nsize,nsize)= zbus(m,m)+zline(i); end % type 3 if (type(i)==3) zsum=zbus(m,m)+zline(i); for k=2:nsize for j=2:nsize zmod(k,j)=-1*zbus(k,m)*zbus(m,j)/zsum; end end zbus=zbus+zmod; end % type 4 if (type(i)==4) zsum=zbus(m,m)+zbus(p,p)-2*zbus(m,p)+zline(i); for k=2:nsize for j=2:nsize zmod(k,j)=-1*(zbus(k,m)-zbus(k,p))*(zbus(m,j)-zbus(p,j))/zsum; end end zbus=zbus+zmod; end display(zbus); end display(zbus);
