AMITY UNIVERSITY
DEPARTMENT OF ELECTRICAL & ELECTRONICS ENGINEERING
MANUAL ON
POWER SYSTEM ANALYSIS THROUGH MATLAB For Mtech (Power Systems) - Ist Year Prepared By:
Ms Shruti Saggar Lecturer ASET
Amity University, Noida
CONTENTS S.No
1.
EXPERIMENTS
Write MATLAB Program to compute performance of Transmission lines. i.e. computation of ABCD parameters for the following models. i)
ii) iii) 2. 3.
Short Transmission Line Medium Line ( Nominal ∏ equivalent ) Lumped parameter ( ∏ equivalent )
Write MATLAB Program to compute sending end voltage and current of transmission line. Write MATLAB Program to compute sending end power, line loss, efficiency and voltage regulation of transmission line.
4.
Write MATLAB Program to plot the voltage profile of a three phase transmission line at the receiving end with varying load conditions.
5. 6.
Write MATLAB program for building Ybus Matrix and use it to find to Ybus matrix for the given power system network. Write Matlab program for power flow studies by Gauss Siedel Method using Ybus.
7. 8.
Write Matlab program for power flow studies by Gauss Siedel Method using Zbus. Write Matlab program for analyzing shunt fault conditions in power system.
EXPERIMENT 1 Write MATLAB Program to compute performance of Transmission lines. i.e. computation of ABCD parameters for the following models. i)
Short Transmission Line Medium Line ( Nominal ∏ equivalent ) Lumped parameter ( ∏ equivalent )
ii) iii)
THEORY : The
transmission line parameters include series resistance and inductance and shunt capacitance. The line models are classified by their length. These classifications are Short line approximation for lines that are less than 80 km long. Medium line approximation for lines whose lengths are between 80 km to 250 km. Long line model for lines that are longer than 250 km.
• • •
The sending and receiving end voltages are denoted by V S and V R respectively. Also the currents I S and I R are entering and leaving the network respectively. The sending end voltage and current are then defined in terms of the ABCD parameters as
A
=
V S
= AV R + BI R
I S
= CV R + DI R
V S
V R
I R
0
=
A is the ratio of sending end voltage to the open circuit receiving B
(
V S
=
I R
end voltage. Similarly,
Ω V R
0
=
.e., B, given in Ohm, is the ratio of sending end voltage and short circuit receiving end current. In a similar way we can also define C
D
=
I S V R
=
mho I R
=
V R
=
0
I S I R
0
Function Used :
The function trmlnper is a generalized MATLAB program for computing the performance of a three phase transmission line. The function computes the ABCD parameter from the line data and provides the option of modeling a) with distributed parameters, b) lumped parameters as i) ∏ equivalent circuit , ii ) nominal ∏ equivalent and iii ) short line Details of the input to the function trmlnper are as follows: r Per phase series resistance of the line Ω per unit length L Per phase inductance of the line in H per unit length g Per Phase shunt conductance of the line in mho per nit length CC Per phase shunt capacitance of the line in F per unit length l length of the line Vr Line to Line voltage in kv at the receiving end Pr three phase receiving end power in MW Pf power factor of the load at the receiving end Pftype type of power factor i.e. lagging or leading. For lagging power factor pftye = 0 During execution, the function trmlnper asks for input in respect of type of line. Type of Line = 1 for lumped parameters when the line is modeled as a ∏ equivalent circuit =2 for medium lines (nominal pi equivalent), = 3 when the line is modeled as a short line. PROGRAM :
Function [A, B, C, D, Vs] = trmlnper(r, L, g, CC, 1, Vr, Pr, pf, pftype); % Programe to compute the performance of transmission line % computing the ABCD constants Vs=zeroes (2, 1); z=r + i*2*pi*50*L; y=g + i*2*pi * 50* CC; gamma = sqrtm (z*y); Zc=sqrtm (z/y); type = input(‘type of line’) % 1 for lumped parameter (pi equivalent), 2 for medium lines (nominal pi equivalent), 3 for short line if type = = 1; Zd = (z * 1*sinh (gamma*1)) / (gamma * 1); Yd=(y*1* tanh (gamma*1/2))/ (gamma* ½); A= (1+ (Yd*Zd/2)); B=Zd; C=Yd*(1+ (Yd*Zd/4)); D=A; else
end if type = = 2; Z=z*1; Y=y*1; A= (1+(Y*Z / 2 )); B=Z; C=Y*(1+(Y*Z/4)); D=A; else end if type = = 3; Z=z*1; A=1; B=Z; C=0; D=A; else end Ir=Pr / (sqrtm (3)*Vr*Pf); Vr=Vr/ sqrtm (3); Phi =acos (pf); If pftye = = 0; phi = -phi; else end
EXPERIMENT 2
Write MATLAB Program to compute sending end voltage and current of transmission line. THEORY: PROGRAM :
Function [A, B, C, D, Vs] = trmlnper(r, L, g, CC, 1, Vr, Pr, pf, pftype); % First compute the ABCD constants (for distributed constants) Vs=zeroes (2, 1); z=r + i*2*pi*50*L; y=g + i*2*pi * 50* CC; gamma = sqrtm (z*y); Zc=sqrtm (z/y); A = cosh (gamma * 1); B= Zc*sinh (gamma*1); C= sinh (gamma * 1)/ Zc; D=A; Ir=Pr / (sqrtm (3)*Vr*Pf); Vr=Vr/ sqrtm (3); Phi =acos (pf); If pftye = = 0; phi = -phi; else end % Computing the Sending end voltage and Current
[Vs] = [A B; C D] * [Vr; Ir] Vs=abs (Vs (1)) deltaVs = angle (Vs (1)) * 180 / pi Is = abs (Vs (2)) deltaIs = angle (Vs (2)) * 180 / pi
EXPERIMENT 3
Write MATLAB Program to compute sending end power, line loss, efficiency and voltage regulation of transmission line. .
Function [A, B, C, D, Vs] = trmlnper(r, L, g, CC, 1, Vr, Pr, pf, pftype); % First compute the ABCD constants (for distributed constants) Vs=zeroes (2, 1); z=r + i*2*pi*50*L; y=g + i*2*pi * 50* CC; gamma = sqrtm (z*y); Zc=sqrtm (z/y); A = cosh (gamma * 1); B= Zc*sinh (gamma*1); C= sinh (gamma * 1)/ Zc; D=A; Ir=Pr / (sqrtm (3)*Vr*Pf); Vr=Vr/ sqrtm (3); Phi =acos (pf); If pftye = = 0; phi = -phi; else end % computing the Sending end voltage and Current [Vs] = [A B; C D] * [Vr; Ir] Vs=abs (Vs (1)) deltaVs = angle (Vs (1)) * 180 / pi Is = abs (Vs (2)) deltaIs = angle (Vs (2)) * 180 / pi
% Computing the Sending end power, line loss, efficiency and voltage regulation
Ps = (Vs (1)* (Vs (2)) ‘) PS=real (Ps)*3 lnloss = PS -Pr Effy = (1- lnloss/PS)*100 Reg = (VS- abs (Vr))/ (abs (Vr)) VS= sqrtm (3)*VS
EXPERIMENT 4
Write MATLAB Program to plot the voltage profile of a three phase transmission line at the receiving end with varying load conditions.
The MATLAB Script file voltage profile is used to compute the voltage profile of the transmission line. Function Used :
Assumptions/Numerical:
600 km long Transmission line Series inductive reactance: XL = 0.35 H/km/phase Shunt Capacitance: Xc = 4.2 µF/km/phase. Assume receiving end voltage equal to 1.0 pu in all cases. Plotting voltage profile at varying loads i. e at 25 %, 50 %, 100 %, 150 % of the SIL. PROGRAM:
% Plotting a Voltage Profile >> L= 0.35 / (2 * pi * 50) L= 0.0011 >> C = 1.3369e-008 >> Zc = 288.6751 >> beta = 2*pi*50*sqrtm*(L*C) Beta= 0.0012 >>Vr=1.0; X=1.0; >>step=600/20;
% computes the inductance of the line
% computes the shunt capacitance of the line %computes the surge impedance of the line % sets the receiving end voltage=1.0 p.u. % X represents the load as a percentage of SIL. % voltage is computed at 20 points. Step fixes %the increment in x which represents the % distance in km from the receiving end.
>>x=600:-step: 0; >> y = (cos (beta*x) + i*sin (beta*x)/x)* Vr; %expression for computing the voltage at %point x from the receiving end >>plot (x, abs(y),’k’)
% plots the distance along the x axis and %absolute value of the voltage at point x km %from the receiving end in black solid line
Hold on % activates the hold for overlay plots >>X=0.25; % fixes the receiving end loss at 25% of SIL. >>y = (cos (beta*x) + i*sin (beta*x)/X)*Vr ;
>>plot(x, abs(y), ‘k-’)
% Plots the distance along the x-axis and %absolute value of the voltage at point x km %from the receiving end in black solid line.
>>X=0.5; >>y = (cos (beta*x) + i*sin (beta*x)/X)*Vr; >>plot(x, abs(y), ‘k+’) % Plots the distance along the x-axis and %absolute value of the voltage at point x km %from the receiving end in black with a + %sign X=0.75; >>y = (cos (beta*x) + i*sin (beta*x)/X)*Vr; >>plot(x, abs(y), ‘k.’) % Plots the distance along the x-axis and %absolute value of the voltage at point x km %from the receiving end in black with a . %point X=1.25; >>y = (cos (beta*x) + i*sin (beta*x)/X)*Vr; >>plot(x, abs(y), ‘k*’) % Plots the distance along the x-axis and %absolute value of the voltage at point x km %from the receiving end in black with a *. %point X=1.5; >>y = (cos (beta*x) + i*sin (beta*x)/X)*Vr; >>plot(x, abs(y), ‘k*’) % Plots the distance along the x-axis and %absolute value of the voltage at point x km %from the receiving end in black with a %diamond >> Xlabel (‘Distance from receiving end in km.’) % labels the x-axis as distance from % receiving end in km. >> Ylabel (‘Sending end voltage in p.u.’)
% labels the y-axis as sending end voltage %in p.u >> Title (‘Voltage profile of a three phase transmission line ’) % provides the title voltage % profile of a three phase transmission line >> hold off
% deactivates the hold for overlay of plots. EXPERIMENT 5
Write MATLAB program for building Ybus Matrix and use it to find to Ybus matrix for the given power system network.
Since Power System networks are composed of linear bilateral components, Ybus is a square symmetric matrix. If elements connected to a bus are not mutually coupled to other elements then the formulation of Ybus matrix of a network is a two step procedure : Theory:
Step 1 : To obtain the diagonal terms, add the admittance of the branches to the ith node. Step 2 : To obtain the off diagonal term , Yij connected between the node I and node j , take the negative of the sum of all the primitive admittance of the branches connected between nodes i and j. Since the matrix is symmetric Yij=Yji. Formulation of Ybus through MATLAB is explained below. Building a Ybus matrix by inspection is very time consuming. Hence using computers for building Ybus matrices is essential. The following is an executable MATLAB programme for building Ybus matrix from line impedance data. Function Used:
The programme is executed through a function file. The programme is saved in a file named admitnat1.m and the output gives the element of Ybus. The various input parameters are as follows. loc: a row vector which stores the start of each row of Ybus elemz: a row vector which stores the primitive impedances of the lines col : a row vector which provides the column number of Ybus nbus : number of buses in the system locs : gives the size of the loc vector. PROGRAM:
Function [Ybus] = admitmat1 (loc, elemz, col, nbus, locs); %Part 1: Program for building Ybus Ybus= zeroes (nbus); locs = locs – 1; nbus = nbus-1; I=1; k=1; for I I = 1: locs; add= loc (I I+1) – loc (I I); for kk=1: add; J=col (k); Ybus (I,I)= Ybus(I,I) + 1/elemz(k); If J = = 0; Disp (‘branch between node & ground’)
k=k+1; else Ybus (J, J) =Ybus (J, J) + 1/elemz(k); Ybus (I, J) =Ybus (I,J) - 1/elemz(k); Ybus (J, I) =Ybus (I, J); k=k+1; end end I=I+1; end Ybus EXPERIMENT 6
Write Matlab programme for power flow studies by Gauss Siedel Method using Ybus THEORY : Load-flow studies are probably the most common
of all power system analysis calculations. They are used in planning studies to determine if and when specific elements will become overloaded. The load-flow problem is as follows: ● At any bus there are four quantities of interest: │V│, θ, P, and Q. ● If any two of these quantities are specified, the other two must not be specified otherwise we end up with more unknowns than equations. ● Because records enable the real and reactive power to be accurately estimated at loads, P and Q are specified quantities at loads, which are called PQ buses. ● Likewise, the real power output of a generator is controlled by the prime mover and the magnitude of the voltage is controlled by the exciter, so and P and │V│ are specified at generators, which are called PV buses. ● This means that │V│ and θ are unknown at each load bus and θ and Q are unknown at each generator bus. ● Since the system losses are unknown until a solution to the load-flow problem has been found, it is necessary to specify one bus that will supply these losses. This is called the slack (or swing , or reference) bus and since P and Q are unknown, │V│ and θ must be specified. Usually, an angle of θ = 0 is used at the slack bus and all other bus angles are expressed with respect to slack. There are two major iterative techniques: The Gauss-Seidel method is based on substituting nodal equations into each other. It is the slower of the two but is the more stable technique. It’s convergence is said to be Monotonic. The iteration process can be visualized for two equations: A ) The Gauss-Seidel Method :
The Power Flow studies using Gauss Siedel Method through MATLAB is explained as below.
Function used:
pfsg
The inputs to the programme functions pfsg are nbus : nlns: ld :
Number of buses in the power system Number of lines Line data. The line data is in the form of a matrix whose dimensions are (nlns x 4 ). The first two columns in the ld matrix are the bus numbers between which the line is connected, the third columns shows the complex impedances of the line and the fourth column provides one half of the line susceptance. The complex impedance and susceptance pertaining line are shown in pu. pfbd: pfbd matrix includes data related to the loads on each bus and the type of buses, that is, slack, load or voltage controlled bus. The dimensions of the matrix are (nbus X 5 ). Column 1 shows the type of bus, 0 for load bus and 2 for a voltage controlled bus. For a slack bus, columns 2 to 4 are zero. For a load bus, columns 2 and 3 indicate real and reactive generated power while columns 4 and 5 show real and reactive load. Lagging reactive load is shown as negative. For a voltage controlled bus , column 2 gives the real power , column 3 shows the specified voltage magnitude and columns 4 and 5 indicate the specified maximum and minimum Q, respectively. accrucy : Tolerance limit Program:
Function [Ybus, VN] = pfsg (nbus, nlns, ld, pfbd, accrucy): %Program for performing power flow studies % Part 1 Ybus = zeroes (nbus) VP (nbus, 1) = 0 VN (nbus, 1) = 0 For k =1: nlns; m=ld(k,1); n=ld(k, 2); Ybus(m,m)=Ybus(m,m) + 1/ld(k,3)+ld(k,4); Ybus(n,n)=Ybus(n,n) + 1/ld(k,3) + ld(k,4); Ybus(m,n)=- 1/ld(k,4); Ybus(n, m) = Ybus(m, n) end Ybus itr =o; Count=0; %Gauss Siedel Iteration Method VP(1) = input (‘slack Bus Voltage ); For k=2: nbus; If pfbd (k, 1) = = 0; VP (k) = VP (1);
else VP (k) = pfbd (k,3); end else while count <=nbus-1; sum = 0; itr = itr+1; for k=2: nbus; if pfbd (k, 1)= = 0 I(k)=(pfbd(k,2)-pfbd(k,4)+i*(pfbd(k,3)-pfbd(k,5)))/(VP(k))’; else sum = 0; for m=1: nbus sum= sum+Ybus (k, m)* VP(m); end q(k) = imag(VP(k)*(sum)’); I(k)= ((pfbd(k,2)+i*q(k))/(VP(k)))’; end sum=0; for m=1:nbus; If m=k sum=sum+Ybus(k,m)*VP(m); else end end VN(k)=(I(k)-sum) / Ybus(k,k); If abs(real(VN(k))-real(VP(k)))<=accrucy; count=count+0.5; else end if pfbd(k,1)= = 0; VP(k)=VN(k); else z=imag(VN(k)); zz=(pfbd(k,3)^2-(z^2))^0.5; VN(k)=zz+i*z; VP(k)= VN(k); end end end fprintf(‘%6.2f\’,itr); fprintf(‘%13.4f\’,real(VP(k)))); fprintf(‘%15.4f\n’,imagl(VP(k)))); end EXPERIMENT 7
Write Matlab program for power flow studies by Gauss Siedel Method using Zbus. Function used:
The MATLAB Program function pfsz computes the bus voltages by the Gauss Iterative method in the Zbus frame of reference. The Zbus matrix for the system employs the MATLAB Program function zeebus. zeebus : The Matlab program function [Zbus]= zeebus(nelemnt,ind)= assembles the
impedance matrix [Zbus] from a partially assembled [Zbus] or from the very beginning. The input parameters are nelemnt for the total number of elements in the network and ind is an indicator. ind =1 if partial [Zbus] is available and ind = 0 if [Zbus] is to be assembled from the beginning. The inputs [nbus, pfbd, accrucy] used for the program function pfsz are similar to the inputs used for the program function pfsg. Function [VP, VN]=pfsz (Zbus, nbus, pfbd, accrucy); % Program for performing power flow studies using Zbus %Part 1 – Output from the function zebus give zbus VP=zero (nbus, 1); VN=zero (nbus, 1); I=zero (nbus, 1) itr=0; count=0; Part 2 : Gauss Iterative Method VP(1)=input (‘slack bus voltage’) For k =2:nbus; VP(k)=VP(1); end while count <=nbus-1; itr = itr+1 for k =2 : nbus; I(k)=(pfbd(k,2)-pfbd(k,4)+i*(pfbd(k,3)-pfbd(k,5)))/(VP(k))’; end for k =2 : nbus; sum = 0; for m =2:nbus; sum =sum + zbus(k-1,m-1)*I(m); end VN(k)=sum+ VP(1) if abs(real(VN(k))-real(VP(k)))
count = count + 0.5; else end VP (k)=VN(k); end fprintf(‘ itr Bus Voltage Voltage \ n’) for k=1:nbus; fprintf(‘%4i \ ‘,itr); fprintf(‘%4i \ ‘,k); fprintf(‘%15.5f\’, real(VP(k))); fprintf(‘%15.5f\’, imag(VP(k))); end end
EXPERIMENT 8 Write Matlab program for analyzing shunt fault conditions in power system.
The MATLAB Program given here is an executable file called faultcal.m. The first line in the program is the function definition line. The variables with in the [ ] define the output variables and the variable within ( ) define the input variables. In the present case, the output variables are the phase and the line voltages, phase currents, and sequence components of voltages and current. Function Used:
The input variables consist of sequence, fault and ground impedances and voltage at the fault point. The input variable type specifies the fault to be analyzed ie that is for SLG, type=1; for LL, type=2; for 2LG, type=3; and for 3Ø, type = 4. function [Ia,Ib,Ic,Va,Vb,Vc,Vab,Vbc,Vca,Iaz,Ia1,Ia2,Vaz,Va1,Va2] =faultcal (type, Vf, Zz, Z1,Z2,Zf ,……Zg) % Program for analyzing shunt fault conditions in power system. a = (-0.5+0.866i); A=[1 1 1; 1 a^2 a;1 a a^2]’ If type = = 1 Iaz = Vf / (Zz +Z1+ Z2 +3 * Zf); Ia1 = Vf / (Zz +Z1+ Z2 +3 * Zf); Ia2 = Vf / (Zz +Z1+ Z2 +3 * Zf); else end If type = = 1 Iaz = 0; Ia1 =Vf / (Z1+ Z2 + Zf); Ia2 =-Ia1 else end If type = = 3 Ia1 = Vf / (Z1+ Zf) + ((Z2+Zf ) * (Zz+Zf+3*Zg)/(Z2+Zz+2*Zf+3*Zg))); Ia2 = - (Zz+Zf+3*Zg)* Ia1 / )/(Z2+Zz+2*Zf+3*Zg); Iaz = -(Z2+Zf)*Ia1)/(Z2+Zz+2*Zf+3*Zg); else end If type = = 4 Ia1 =Vf/(Z1+Zf); Ia2 = 0; Iaz = 0; else end Vaz=Zz*Iaz; Va1=Vf-Z1*Ia1; Va2=-Z2*Ia2; I= [ Iaz; Ia1 ; Ia2 ]; V= [Vaz; Va1; Va2 ]; I I = A* 1; Ia= I I (1)
Ib= I I(2) Ic= I I ( 3 ); W=A* V ; Va=VV(1); Vb=VV(2); Vc=VV(3); Vab=Va-Vb; Vbc=Vb-Vc; Vca=Vc-Va
