Experiment Introduction to powersystems学生姓名：学号：专业班级：实验名称：电力系统导论（双语）CONTENTS1、EXPERIMENT 1 (1)B US A DMITTANCE M A TRIX ...................................................................................................................... 1-82、EXPERIMENT 2 (8)B US I MPEDANCE M A TRIX....................................................................................................................... 8-213、EXPERIMENT 3 (21)G AUSS-S EIDEL AND NEWTON METHOD................................................................................................. 21-244、PERSENAL SUMMARY (24)1 Experiment 1Bus Admittance Matrix1. Objective• To write a simple program in MATLAB® for the algorithm of bus admittance matrix. 2. DiscussionBus Admittance matrixThe node-voltage equation of an n-bus power system can be written in matrix form as(1)Or（2）Where Ibus is the column vector of the injected bus currents. Vbusis the column vector of busvoltages measured from the reference node. Y busis known as the bus admittance matrix .Calculation method:1. The diagonal element of each node is the sum of admittances connected to it. It is knownas the self-admittance or driving point admittance , i.e.,（3）2. The off-diagonal element is equal to the negative of the admittance between the nodes. Itis known as the mutual admittance or transfer admittance , i.e.(4)Result:Base on t he equations (3) and (4), the algorithm of bus admittance matrix can be used to build the bus admi ttance matrix (Y bus). 3. Mathmatics modelThe mathematics models of power system elements are two kinds, one is transmission line or cable and another is transformer, and they are equivalent into pi type equivalent circuit, and2described as following: 1) Transmission line or cableAccording to the Pi sectional circuit of transmission line or cable, we can get the elements of nodal Admittance matrix value, as the followingFigure: The pi section equivalent circuit of transmission line or cable()ij ijio ii ii jX Rg Y Y +++=1(5) ()ij ijjo jj jj jX Rg Y Y +++=1(6)()ij ijij ij jX RY Y +-=1(7)ijji Y Y =(8)2) transformerAs show as the transformer circuit , we have k I Z V V jT j i -=(9)j i KI I -=(10)According to equations (10) and (9) , we have()()()i j Tj T j T i j T j T i j T j V V Z K V Z K V Z K V V Z K V Z KV V Z I -+-=--+=-=111(11)()()i Tj i T i T j T i T j i T j i V Z K V V Z K V Z K V Z K V Z K V KV Z K KI I --+=-=-=-=22()j i Ti T V V Z KV Z K K -+-=2(12)j3Figure: the Pi section equivalent circuit of TransformerAssume ijij T T jX R Z Y +==11 ()ij ijii ii jX R KY Y ++=2(13) ()ij ijjj jj jX RY Y ++=1(14) ()ij ijij ij jX RKY Y +-=(15)ijji Y Y =(16)4. System RequirementComputer with MATLAB® 6 or above installed. 5. Procedure1.0 Launch the MATLAB program.2.0 Go to FILE NEW M-file.3.0 Write a function Y = The_Node_Admittance_Matrix(TopoStructureAndBranchPara)for the formation of the bus admittance matrix.4.0 TopoStructureAndBranchPara is the transmission line, cable and transformer inputdata and contains five columns parameters. The first two columns are the line bus numbers and the remaining columns contain the line resistance and reactance in per-unit and transformer tap ratio or capacitor of transmission line.5.0 The function should return the bus admittance matrix.6. ExercisesUse the written function, Y = The_Node_Admittance_Matrix (TopoStructureAndBranchPara) to obtain the Y bus of the following power system network:Q1. You are required to write the Ybus topological structure and parameter into a text file.(Hint: use the matlab text compiler to write down the table 1 data, using the comma to separate the parameters, and save it use the name of 4_Power_System_Data.dbf)Q2. You are required to write out the program flow figure of forming a nodal admittance matrix.Hint. You are required to compile a program to form the Y bus Matrix, the following program isa reference program to you.Figure : One-line diagram of power system47.The flow chartFigure : The flow chart of Forming Nodal Admittance Matrix5程序是：%function OutPut=The_Node_Admittance_Matrix(handles)%is a subroutine of PowerSystemCalculationfunction OutPut=The_Node_Admittance_Matrix(handles)%the following program is open a data file and get the Number of% Node and Branch data to form a nodal addmittance matrix%the following code is open a file and read the data of power system network[fname,pname] = uigetfile('*.dbf','Select the network parametre data-file'); TopoStructureAndBranchPara= csvread(fname);[NumberOfBranch,NumberOfPara]=size(TopoStructureAndBranchPara);Temporary1=max(TopoStructureAndBranchPara(:,1));Temporary2=max(TopoStructureAndBranchPara(:,2));if Temporary1 > Temporary2NumberOfNode=Temporary1;elseNumberOfNode=Temporary2;end%The following program is to form the Nodal Admittance Matrix% and the Topologic structure and Branch Parametres are arranged% I,J,R,X,C/K, and pay attention to the inpedence of transformer is in the% side of Node J and the ratio of transformer 1:K is in the side of Node Ifor CircleNumber1=1:NumberOfBranchfor CircleNumber2=1:NumberOfBranchNodalAdmittanceMatrix(CircleNumber1,CircleNumber2)=0;endendfor CircleNumber=1:NumberOfBranchif TopoStructureAndBranchPara(CircleNumber,5) > 0.85NodalAdmittanceMatrix(TopoStructureAndBranchPara(TopoStructureAndBranchPara(CircleNum ber,1),TopoStructureAndBranchPara(CircleNumber,1)))=...NodalAdmittanceMatrix(TopoStructureAndBranchPara(TopoStructureAndBranchPara(CircleNum ber,1),TopoStructureAndBranchPara(CircleNumber,1)))+...TopoStructureAndBranchPara(CircleNumber,5)^2/...(TopoStructureAndBranchPara(CircleNumber,3)+...j*TopoStructureAndBranchPara(CircleNumber,4)) ;NodalAdmittanceMatrix(TopoStructureAndBranchPara(CircleNumber,2),TopoStructureAndBranc hPara(CircleNumber,2))=...NodalAdmittanceMatrix(TopoStructureAndBranchPara(CircleNumber,2),TopoStructureAndBranc hPara(CircleNumber,2))+...61/((TopoStructureAndBranchPara(CircleNumber,3)+j*TopoStructureAndBranchPara(CircleNumb er,4)));NodalAdmittanceMatrix(TopoStructureAndBranchPara(CircleNumber,1),TopoStructureAndBranc hPara(CircleNumber,2))=...NodalAdmittanceMatrix(TopoStructureAndBranchPara(CircleNumber,1),TopoStructureAndBranc hPara(CircleNumber,2))...-TopoStructureAndBranchPara(CircleNumber,5)/...((TopoStructureAndBranchPara(CircleNumber,3)+j*TopoStructureAndBranchPara(CircleNumber, 4)));NodalAdmittanceMatrix(TopoStructureAndBranchPara(CircleNumber,2),TopoStructureAndBranc hPara(CircleNumber,1))=...NodalAdmittanceMatrix(TopoStructureAndBranchPara(CircleNumber,1),TopoStructureAndBranc hPara(CircleNumber,2));elseNodalAdmittanceMatrix(TopoStructureAndBranchPara(CircleNumber,1),TopoStructureAndBranc hPara(CircleNumber,1))=...NodalAdmittanceMatrix(TopoStructureAndBranchPara(CircleNumber,1),TopoStructureAndBranc hPara(CircleNumber,1))+...+1/(TopoStructureAndBranchPara(CircleNumber,3)+...j*TopoStructureAndBranchPara(CircleNumber,4))+j*TopoStructureAndBranchPara(CircleNumbe r,5);NodalAdmittanceMatrix(TopoStructureAndBranchPara(CircleNumber,2),TopoStructureAndBranc hPara(CircleNumber,2))=...NodalAdmittanceMatrix(TopoStructureAndBranchPara(CircleNumber,2),TopoStructureAndBranc hPara(CircleNumber,2))+...+1/(TopoStructureAndBranchPara(CircleNumber,3)+...j*TopoStructureAndBranchPara(CircleNumber,4))+j*TopoStructureAndBranchPara(CircleNumbe r,5)NodalAdmittanceMatrix(TopoStructureAndBranchPara(CircleNumber,1),TopoStructureAndBranc hPara( CircleNumber,2))=...NodalAdmittanceMatrix(TopoStructureAndBranchPara(CircleNumber,1),TopoStructureAndBranc7hPara( CircleNumber,2))...-1/(TopoStructureAndBranchPara(CircleNumber,3)+...j*TopoStructureAndBranchPara(CircleNumber,4));NodalAdmittanceMatrix(TopoStructureAndBranchPara( CircleNumber,2),TopoStructureAndBran chPara(CircleNumber,1))=...NodalAdmittanceMatrix(TopoStructureAndBranchPara(CircleNumber,1),TopoStructureAndBranc hPara( CircleNumber,2));endend运行结果:NodalAdmittanceMatrix =1.0421 - 8.2429i -0.5882 +2.3529i 0 +3.6667i -0.4539 + 1.8911i-0.5882 + 2.3529i 1.0690 - 4.7274i 0 00 + 3.6667i 0 0 - 3.3333i 0-0.4539 + 1.8911i 0 0 0.9346 - 4.2616i Experiment 2Bus Impedance Matrix1. Objective• To write a simple program in MATLAB® for the algorithm of bus impedance matrix.2. DiscussionBus Impedance matrixThe node-current equation of an n-bus power system can be written in matrix form as8(1) Or（2）Where Ibus is the vector of the injected bus currents. Vbusis the vector of bus voltagesmeasured from the reference node. Zbusis known as the bus impedance matrix.Z bus CALCULATIONThere are two major methods:1. Inversion of Y bus2. Z bus building algorithm3. LU factorization.Matrix inversion for very large-scale networks is very slow and avoided. On the other hand, the Z bus building algorithm is the acceptable method for large networks. It is much faster and easier to update.Z bus Building AlgorithmThis is a very systematic way of forming Z bus and allow for easy implementation on the digital computer.Advantages1. Easy to implement2. Avoids matrix inversion3. It is easy to handle Z bus modifications due to line switchings and changes in network topologySolution：NotationThe addition of a branch of impedance z b falls into 4 categories as follows. CASE 1:Adding z b from a new bus p to the reference bus 0Here, ●new Bus Z will be ()()11++n n matrix● Injected current p I will not change the bus voltages of the original network ●p b p I z V =Then,CASE 2:Adding z b from a new bus p to an existing bus kHere, ●new Bus Z will be ()()11++n n matrix● Injected current p I will change the bus voltages of the original network ● Injected current at bus k becomes p k I I + ●p b k p I z V V +=Then,Summary of Case 2CASE 3: Adding z b from an existing bus k to the reference bus 0Here, ●new Bus Z will be n n ⨯ matrix● This is a special case of CASE 2 except that 0=p VThen,Since Vp = 0, the last row can be eliminated by Kron reduction (one kind of mathematics method) to yield an n×n matrix asSummary of Case 3 • Proceed as in CASE 2 • Eliminate the last rowCASE 4: Adding z b from two existing buses j and kHere,●new Bus Z will be n n ⨯ matrix● Branch current p I will change the bus voltages of the original network ● Injected current at bus j becomes p j I I + ● Injected current at bus k becomes p k I I - ●p b k j p b j k I z V V I z V V +-=→=-0● will change the bus voltages of the original network ● Injected current at bus k becomes p k I I +p b k p I z V V +=Then,The last row can be eliminated by Kron reduction to yield an n×n matrix asSummary of Case 4 1. Proceed as follows2. Eliminate the last rowPROCEDURE FOR BUILDING ALGORITHMSteps1. Prepare a transmission line list●First line considered must have a tie to the reference●All subsequent lines in the list must have at least one bus already seen2．Algorithm●Determine which of the cases (1-4) apply and proceed accordingly.●Read a line. If the last line has been processed, then stop else consideranother lineComments on the algorithm1. Loops must be closed as early as possible2. The axes of the final matrix may not correspond with the buses in numerical order3. The line list order affects the efficiency of calculation4. If there is no tie to the reference, place a fictitious ground tie at one of the buses.5. Line outage is modeled by using –z b instead of z b3. System RequirementComputer with MATLAB® 6 or above installed.4. Procedure1.0 Launch the MATLAB program.2.0 Go to FILE NEW M-file.3.0 Write a function Z = znbus (z) for the formation of the bus impedance matrix.4.0 z is the line input and contains three columns. The first two columns are the line busnumbers and the remaining columns contain the line resistance in per-unit.5.0 The function should return the bus impedance matrix.5. Mathmatics modelThe mathmatics models of power system elements are two kinds, one is transmission line and another is transformation, and they are equivalented into pi type equivalent ciucuit,and dricribled as following:6. ExercisesUse the written function, Z = znbus(z) to obtain the Ybus of the following power system network:Example 1Figure 3: One-line diagram of power systemQ2. You are required to write the Z bus into a text file. (Hint: use the matlab text compiler) Example 2For the system shown, form Zbus matrix using the building algorithmjij1:KSolutionA line listApply Kron reduction to eliminate the last rowHint. You are required to compile a program to form the Z bus Matrix.the following programis a reference program to you.The program is:%function OutPut=The_Node_impedance_Matrix(handles)%is a subroutine of PowerSystemCalculationfunction OutPut=The_Node_impedance_Matrix(handles)%the following program is open a data file and get the Number of% Node and Branch data to form a nodal impedance matrix%the following code is open a file and read the data of power system network[fname,pname] = uigetfile('*.dbf','Select the network parametre data-file');Topo_Structure_And_Branch_Para= csvread(fname);%get the electric power system the number of branch and the parametre of% elements[NumberOfBranch,NumberOfPara]=size(Topo_Structure_And_Branch_Para);%Temporary1---temporary variable 1%Temporary2---temporary variable 2Temporary1=max(Topo_Structure_And_Branch_Para(:,1));Temporary2=max(Topo_Structure_And_Branch_Para(:,2));if Temporary1 > Temporary2NumberOfNode=Temporary1;elseNumberOfNode=Temporary2;end% The following program is to form the Nodal impedance Matrix% and the Topologic structure and Branch Parametres are arranged% I,J,R,X,C/K, and pay attention to the inpedence of transformer is in the% side of Node J and the ratio of transformer 1:K is in the side of Node %% set the initial value of Nodal Admittance Matrix to zerofor CircleNumber1=1:NumberOfNodefor CircleNumber2=1:NumberOfNodeNodal_impedance_Matrix(CircleNumber1,CircleNumber2)=0;endendfor CircleNumber=1:NumberOfBranchif Topo_Structure_And_Branch_Para(CircleNumber,5) > 0.85Nodal_impedance_Matrix(Topo_Structure_And_Branch_Para(Topo_Structure_And_Branch_Para (CircleNumber,1),Topo_Structure_And_Branch_Para(CircleNumber,1)))=...Nodal_impedance_Matrix(Topo_Structure_And_Branch_Para(Topo_Structure_And_Branch_Para (CircleNumber,1),Topo_Structure_And_Branch_Para(CircleNumber,1)))+Topo_Structure_And_B ranch_Para(CircleNumber,5)^2/(Topo_Structure_And_Branch_Para(CircleNumber,3)+...j*Topo_Structure_And_Branch_Para(CircleNumber,4)) ;Nodal_impedance_Matrix(Topo_Structure_And_Branch_Para(CircleNumber,2),Topo_Structure_ And_Branch_Para(CircleNumber,2))=...Nodal_impedance_Matrix(Topo_Structure_And_Branch_Para(CircleNumber,2),Topo_Structure_ And_Branch_Para(CircleNumber,2))+...1/((Topo_Structure_And_Branch_Para(CircleNumber,3)+j*Topo_Structure_And_Branch_Para(Ci rcleNumber,4)));Nodal_impedance_Matrix(Topo_Structure_And_Branch_Para(CircleNumber,1),Topo_Structure_ And_Branch_Para(CircleNumber,2))=...Nodal_impedance_Matrix(Topo_Structure_And_Branch_Para(CircleNumber,1),Topo_Structure_ And_Branch_Para(CircleNumber,2))...-Topo_Structure_And_Branch_Para(CircleNumber,5)/...((Topo_Structure_And_Branch_Para(CircleNumber,3)+j*Topo_Structure_And_Branch_Para(Circ leNumber,4)));Nodal_impedance_Matrix(Topo_Structure_And_Branch_Para(CircleNumber,2),Topo_Structure_ And_Branch_Para(CircleNumber,1))=...Nodal_impedance_Matrix(Topo_Structure_And_Branch_Para(CircleNumber,1),Topo_Structure_ And_Branch_Para(CircleNumber,2));elseNodal_impedance_Matrix(Topo_Structure_And_Branch_Para(CircleNumber,1),Topo_Structure_ And_Branch_Para(CircleNumber,1))=...Nodal_impedance_Matrix(Topo_Structure_And_Branch_Para(CircleNumber,1),Topo_Structure_ And_Branch_Para(CircleNumber,1))+...+1/(Topo_Structure_And_Branch_Para(CircleNumber,3)+...j*Topo_Structure_And_Branch_Para(CircleNumber,4))+j*Topo_Structure_And_Branch_Para(Cir cleNumber,5);Nodal_impedance_Matrix(Topo_Structure_And_Branch_Para(CircleNumber,2),Topo_Structure_ And_Branch_Para(CircleNumber,2))=...Nodal_impedance_Matrix(Topo_Structure_And_Branch_Para(CircleNumber,2),Topo_Structure_ And_Branch_Para(CircleNumber,2))+...+1/(Topo_Structure_And_Branch_Para(CircleNumber,3)+...j*Topo_Structure_And_Branch_Para(CircleNumber,4))+j*Topo_Structure_And_Branch_Para(Cir cleNumber,5)Nodal_impedance_Matrix(Topo_Structure_And_Branch_Para(CircleNumber,1),Topo_Structure_ And_Branch_Para( CircleNumber,2))=...Nodal_impedance_Matrix(Topo_Structure_And_Branch_Para(CircleNumber,1),Topo_Structure_ And_Branch_Para( CircleNumber,2))...-1/(Topo_Structure_And_Branch_Para(CircleNumber,3)+...j*Topo_Structure_And_Branch_Para(CircleNumber,4));Nodal_impedance_Matrix(Topo_Structure_And_Branch_Para( CircleNumber,2),Topo_Structure_ And_Branch_Para(CircleNumber,1))=...Nodal_impedance_Matrix(Topo_Structure_And_Branch_Para(CircleNumber,1),Topo_Structure_ And_Branch_Para( CircleNumber,2));endendformat shortNodal_impedance_Matrix*inv(Nodal_impedance_Matrix)运行结果：Nodal_impedance_Matrix =1.0421e+000 -8.2429e+000i -5.8824e-001 +2.3529e+000i 0 +3.6667e+000i 0-5.8824e-001 +2.3529e+000i 5.8824e-001 -2.3377e+000i 0 00 +3.6667e+000i 0 0 -3.3333e+000i 00 0 0 4.5386e-001 -1.8719e+000iNodal_impedance_Matrix =1.0421e+000 -8.2429e+000i -5.8824e-001 +2.3529e+000i 0 +3.6667e+000i -4.5386e-001 +1.8911e+000i-5.8824e-001 +2.3529e+000i 1.0690e+000 -4.7274e+000i 0 00 +3.6667e+000i 0 0 -3.3333e+000i 0-4.5386e-001 +1.8911e+000i 0 0 9.3463e-001 -4.2616e+000ians =1.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i-0.0000 - 0.0000i 1.0000 - 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i-0.0000 - 0.0000i -0.0000 - 0.0000i 1.0000 - 0.0000i -0.00000 - 0.0000i 0 + 0.0000i 0.0000 - 0.0000i 1.0000 + 0.0000i以上就是对阻抗矩阵的验证，其和其逆相乘为单位对角矩阵3、Experiment 3Gauss-Seidel Method1. Objective• To write a simple program in MATLAB® for the algorithm to solution of nonlinear algebraic equations;• Known as the method of successive displacements.2. DiscussionThe most common methods for solving nonlinear algebraic equations are Gauss-Seidel,Newtow-Rahpson, and quasi-Newton-Raphson methods. We start with one dimensional equations and then generalize to n-dimensional equations.3. Mathmatics modelConsider the nonlinear equation 0)(=x f .The equation is broken into two parts thus:)(x g x =. We assume )0(x is an initial "guess" of the solution, then "refine" the solution using:)()0()1(x g x =This process is repeated thus)()1()2(x g x =and on the th n iteration we have: )()1()(-=n n x g x. If this process is convergent, then the successive solutions approach a value which is declared as the solution. Thus if at some step 1+k we have:ε≤-+)()1(k k x xwhere e ε is the desired "accuracy", then we claim the solution has been found to the accuracy specified.4. System RequirementComputer with MATLAB® 6 or above installed.5. Procedure1.0 Launch the MATLAB program.2.0 Go to FILE NEW M-file.3.0 Write a function program of Gauss Seidel Method.6. Exercises Example: Using the Gauss-Seidel method to obtain the roots of the equation:0496)(23=-+-=x x x x f First the equation is expressed in a different form thus())(469123x g x x x =---= And the iteration can proceed. Take a good look at the shape of the iterations! Below is the program showing the process graphically (later showing how to do it iteratively).7.The flow chart of Gauss Seidel method (Omitted)8.Reference Program and result.程序是：clear allclcx0=0.5;n=0;while (abs(x0^3-6*x0^2+9*x0-4)>0.00001)x0=-(x0^3-6*x0^2-4)/9;y=x0;n=n+1;end结果是：n=1627 y=x0=0.99818clear allclcx0=2.5;n=0;while (abs(x0^3-6*x0^2+9*x0-4)>0.00001)x0=-(x0^3-6*x0^2-4)/9;y=x0;n=n+1;end结果是：n=7 y=x0=4仿照高斯--赛德尔法，我们可以写出简单的牛顿法的程序,如下：牛顿法解方程x0=0.5;n=0;while (abs((x0^3-6*x0^2+9*x0-4)/(3*x0^2-12*x0+9))>0.00001)dx0=-(x0^3-6*x0^2+9*x0-4)/(3*x0^2-12*x0+9);x0=x0+dx0;n=n+1;end结果是：dx0=1.7684e-005 n=15 x0=0.99998 y= -0.875x0=0.5;n=0;while (abs(x0^3-6*x0^2+9*x0-4)>0.00001)dx0=-(x0^3-6*x0^2+9*x0-4)/(3*x0^2-12*x0+9);x0=x0+dx0;n=n+1;end结果是: dx0=0.0011305 n=9 x0=0.99887 y= -0.875x0=3.5;n=0;while (abs(x0^3-6*x0^2+9*x0-4)>0.00001)dx0=-(x0^3-6*x0^2+9*x0-4)/(3*x0^2-12*x0+9);x0=x0+dx0;n=n+1;end结果是：dx0= -2.5283e-006 n=5 x0=4 y= -0.875实验总结：在做实验前,我以为不会难做,就像以前做物理实验一样,做完实验,然后两下子就将实验报告做完.直到做完测试实验时,我才知道其实并不容易做,但学到的知识与难度成正比,使我受益匪浅.在做实验前,一定要将实验指导书上的知识吃透,因为这是做实验的基础,否则,在老师讲解时就会听不懂,这将使你在做实验时的难度加大,浪费做实验的宝贵时间。