1、实验报告06级昆明理工大学 电力工程 学院学生实验报告实验课程名称： 电力系统计算机辅助分析上机实验 开课实验室： 七号楼计算机中心机房 2009年10月19日年级、专业、班电自061学号姓名成绩实验项目 名 称电力系统计算中常用的数值算法指 导 教 师宋琪教 师 评 语 教师签名：宋琪年 月 日一、上机实验目的1、了解、掌握MATLAB软件的基本运用2、掌握高斯消去法、因子表求解线性方程组的基本原理3、读懂课本中相应算法的实际例程4、掌握MATLAB针对上述算法提供的函数二、算法功能与原理 三、程序流程图图17.1形成因子表的框图四、过程记录例：17.1 求以下系数矩阵的因子表例 17.2
2、在例17.1因子表的基础上进行前代输入数据为：请输入矩阵A=1 5 6 0;5 2 7 0;6 7 3 8 ;0 0 8 4输入常数项矩阵B=1 ；2 ；3 ；4结果：利用因子表对常数项进行前代的结果为：B=1.0000 0.1304 0 0.3846例17.3 求例17.2线性方程组的解。输入数据为：请输入矩阵A= 1 5 6 0;5 2 7 0;6 7 3 8;0 0 8 4请输入常数项矩阵B=1； 2 ；3 ；4五、源码1.因子表的形成过程clc;clear all;A=1 5 6 0;5 2 7 0;6 7 3 8;0 0 8 4;n,m=size(A);for i=1:n %按行进行
3、循环 A(i,i)=1./A(i,i); %用对角元位置存储归一化该行所用的系数 for j=i+1:n %用归一化系数乘以本行对角元右侧的元素，并存于原位 A(i,j)=A(i,j)*A(i,i);end for k=i+1:n %用被归一化的行，消去下面的剩余行 for j=i+1:n %实际被处理的元素，在被归一化的对角元的右下角 A(k,j)=A(k,j)-A(k,i)*A(i,j); end endenddisp(矩阵的因子表为：);disp(A)运行结果：矩阵的因子表为： 1.0000 5.0000 6.0000 0 5.0000 -0.0435 1.0000 0 6.0000 -
4、23.0000 -0.1000 -0.8000 0 0 8.0000 0.09622.利用因子表的前代过程clc;clear all;A=1 5 6 0;5 2 7 0;6 7 3 8;0 0 8 4;n,m=size(A);for i=1:n %按行进行循环 A(i,i)=1./A(i,i); %用对角元位置存储归一化该行所用的系数 for j=i+1:n %用归一化系数乘以本行对角元右侧的元素，并存于原位 A(i,j)=A(i,j)*A(i,i);end for k=i+1:n %用被归一化的行，消去下面的剩余行 for j=i+1:n %实际被处理的元素，在被归一化的对角元的右下角 A(
5、k,j)=A(k,j)-A(k,i)*A(i,j); end endenddisp(矩阵的因子表为：);disp(A) B=1 2 3 4;for i=1:4 B(i)=B(i)*A(i,i); %用因子表对角元进行归一化for j=i+1:4 B(j)=B(j)-A(j,i)*B(i);%用因子表对角元下面的元素作为消去系数endenddisp(B); 运行结果：矩阵的因子表为： 1.0000 5.0000 6.0000 0 5.0000 -0.0435 1.0000 0 6.0000 - 23.0000 -0.1000 -0.8000 0 0 8.0000 0.0962 1.0000 0.
6、1304 00.38463.利用因子表的回代过程clc;clear all;A=1 5 6 0;5 2 7 0;6 7 3 8;0 0 8 4;n,m=size(A);for i=1:n %按行进行循环 A(i,i)=1./A(i,i); %用对角元位置存储归一化该行所用的系数 for j=i+1:n %用归一化系数乘以本行对角元右侧的元素，并存于原位 A(i,j)=A(i,j)*A(i,i);end for k=i+1:n %用被归一化的行，消去下面的剩余行 for j=i+1:n %实际被处理的元素，在被归一化的对角元的右下角 A(k,j)=A(k,j)-A(k,i)*A(i,j); en
7、d endenddisp(矩阵的因子表为：);disp(A) B=1 2 3 4;for i=1:4 B(i)=B(i)*A(i,i); %用因子表对角元进行归一化for j=i+1:4 B(j)=B(j)-A(j,i)*B(i);%用因子表对角元下面的元素作为消去系数endenddisp(B); m,n=size(A);for i=n-1:-1:1 %从倒数第二行开始回代，直到第一行 for j=i+1:1:n %取用因子表中，每行对角元右边的元素 B(i)=B(i)-A(i,j)*B(j); endenddisp(B); 运行结果：矩阵的因子表为： 1.0000 5.0000 6.0000
8、 0 5.0000 -0.0435 1.0000 0 6.0000 -23.0000 -0.1000 -0.8000 0 0 8.0000 0.0962 1.0000 0.1304 0 0.3846 0.0401 -0.1773 0.30770.3846六、讨论昆明理工大学 电力工程 学院学生实验报告实验课程名称： 电力系统计算机辅助分析上机实验 开课实验室： 七号楼计算机中心机房 2009年10月26日年级、专业、班电自061学号姓名成绩实验项目 名 称电力网络的数学模型指 导 教 师宋琪教 师 评 语 教师签名：宋琪年 月 日一、上机实验目的1、掌握节点导纳矩阵和节点阻抗矩阵的生成算法2、
9、掌握电力系统计算所需原始数据的准备过程3、读懂课本中相应算法的实际例程二、算法功能与原理 三、程序流程图四、过程记录例 17.6 用节点导纳矩阵的程序求图17.10所示的节点网络的导纳矩阵。解：输入数据：请输入节点数:n=5请输入支路数:nl=5请输入由支路参数形成的矩阵:B=1 2 0.03i 0 1.05 0; 2 3 0.08+0.3i 0.5i 1 0; 3 4 0.015i 0 1.05 1; 2 5 0.1+0.35i 0 1 0; 3 5 0.04+0.25i 0.5i 1 0请输入由节点号及其对地阻抗形成的矩阵:X=1 0;2 0;3 0;4 0;5 0例 17.7 形成如图1
10、7.20所示网络的阻抗阵。解：输入数据如下: 请输入节点数:n=3 请输入支路数:nl=6请输入由支路参数形成的矩阵: B=0 1 2i 0 1 0;0 2 4i 0 1 0;1 2 2i 0 1 0; 0 3 20i 0 1 0;2 3 8i 0 1 0; 1 3 5i 0 1 0 五、源码4.节点导纳矩阵的形成过程clc;clear all;n=5;%节点数nl=5;%支路数B=1 2 0.03i 0 1.05 0;2 3 0.08+0.3i 0.5i 1 0;4 3 0.015i 0 1.05 0;2 5 0.1+0.35i 0 1 0;3 5 0.04+0.25i 0.5i 1 0;d
11、isp(支路参数阵B=);disp(B);X=1 0;2 0;3 0;4 0;5 0;%节点并联支路参数Y=zeros(n);%确定导纳阵阶数，并用进行初始化%对节点上的并联支路进行处理for i=1:n if X(i,2)=0; p=X(i,1); %找到有接地支路的节点号 Y(p,p)=1./X(i,2); %节点接地支路参数直接计入自导纳 end endfor i=1:nl %循环处理原始数据中的每一条支路 % p,q分别代表导纳矩阵中元素的下标 if B(i,6)=0 %变压器参数要求归算至始端，注意第3、支路与1、支路的区别 p=B(i,1); q=B(i,2); else p=B(
12、i,2); q=B(i,1); end %非对角元处理 Y(p,q)=Y(p,q)-1./(B(i,3)*B(i,5); Y(q,p)=Y(p,q); %对角元分为两种情况，p是首端，q是末端%因为变压器阻抗归算至首端,所以从末端看到的变压器阻抗要经过归算 Y(q,q)=Y(q,q)+1./(B(i,3)*B(i,5)2)+B(i,4)./2; Y(p,p)=Y(p,p)+1./B(i,3)+B(i,4)./2;enddisp(导纳矩阵Y=:);disp(Y) 运行结果：支路参数阵B=1.0000 2.0000 0 + 0.0300i 0 1.0500 0 2.0000 3.0000 0.08
13、00 + 0.3000i 0 + 0.5000i 1.0000 0 4.0000 3.0000 0 + 0.0150i 0 1.0500 0 2.0000 5.0000 0.1000 + 0.3500i 0 1.0000 0 3.0000 5.0000 0.0400 + 0.2500i 0 + 0.5000i 1.0000 0 导纳矩阵Y=: 0 -33.3333i 0 +31.7460i 0 0 0 0 +31.7460i 1.5846 -35.7379i -0.8299 + 3.1120i 0 -0.7547 + 2.6415i 0 -0.8299 + 3.1120i 1.4539 -66
14、.9808i 0 +63.4921i -0.6240 + 3.9002i 0 0 0 +63.4921i 0 -66.6667i 0 0 -0.7547 + 2.6415i -0.6240 + 3.9002i 0 1.3787 - 6.2917i5.节电阻抗矩阵的形成clc;clear;n=3;nl=6;B=0 1 2i 0 1 0;0 2 4i 0 1 0; 1 2 2i 0 1 0; 0 3 20i 0 1 0; 2 3 8i 0 1 0; 1 3 5i 0 1 0;m=0;%用于存储最大节点号Z=zeros(n);for k1=1:nl p=B(k1,1); q=B(k1,2); if
15、B(k1,6)=0 %变压器支路首端是低压侧 l=1./B(k1,5); else l=B(k1,5);%变压器支路首端是高压侧 end if p=0 if qm %追加接地树支 Z(q,q)=B(k1,3); m=m+1; else %追加接地连支 for k=1:m, Z(k,m+1)=-Z(k,q); Z(m+1,k)=-Z(q,k); end Z(m+1,m+1)=Z(q,q)+B(k1,3); for l1=1:m for k=1:m Z(l1,k)=Z(l1,k)-Z(l1,m+1)*Z(m+1,k)./Z(m+1,m+1); end Z(l1,m+1)=0; end for k=
16、1:m+1 Z(m+1,k)=0; end end else if qm %追加不接地树支 for k=1:m Z(k,q)=Z(k,p)*l; Z(q,k)=Z(p,k)*l; end Z(q,q)=l2*Z(p,p)+l2*B(k1,3); m=m+1; else for k=1:m %追加不接地连支 Z(k,m+1)=l*Z(k,p)-Z(k,q); Z(m+1,k)=l*Z(p,k)-Z(q,k); end Z(m+1,m+1)=l2*Z(p,p)+Z(q,q)-2*l*Z(p,q)+l2*B(k1,3); for l1=1:m for k=1:m Z(l1,k)=Z(l1,k)-Z(
17、l1,m+1)*Z(m+1,k)./Z(m+1,m+1); end Z(l1,m+1)=0; end for k=1:m+1 Z(m+1,k)=0; end end endendZ=Z(1:n,1:n);disp(阻抗矩阵Z=);disp(Z) 运行结果：阻抗矩阵Z= 0 + 1.4124i 0 + 0.9605i 0 + 1.0734i 0 + 0.9605i 0 + 1.8531i 0 + 1.1299i 0 + 1.0734i 0 + 1.1299i 0 + 3.6158i六、讨论昆明理工大学 电力工程 学院学生实验报告实验课程名称： 电力系统计算机辅助分析上机实验 开课实验室： 七号楼
18、计算机中心机房 2009年11月2日年级、专业、班电自061学号姓名成绩实验项目 名 称电力系统故障的计算机算法指 导 教 师宋琪教 师 评 语 教师签名：宋琪年 月 日一、上机实验目的1、掌握电力系统对称故障计算的基本原理及算法2、掌握电力系统简单不对称故障计算的基本原理及算法3、读懂课本中相应算法的实际例程二、算法功能与原理 三、程序流程图四、过程记录例18.2 图18.6电力系统，负荷全部略去，简化后的各电抗标么值注于等值网络中。试计算f点三相短路时的短路电流及网络中的电流分布。例 18.3 如图18.5所示的网络，母线3发生三相直接短路,试作下列计算: 故障后0秒、2秒、4秒时母线3的
19、短路电流; 故障后2秒、4秒时电源与故障点间的支路电流。各元件的参数如下:发电机: ：100MVA; ：200MVA。额定电压均为10.5KV，次暂态电抗 均为0.2。变压器： ：100MVA; ：200MVA。变比均为10.5/115KV，短路电压百分数均为10。电力线路：三条电力线路 参数均为115KV,60Km ，电抗 =0.44/Km ，电容 负荷： ，50MW，cos=0.985; ， 100MW ， cos=1 。例 18.4 在例题18.1的系统中，已知下列情况：二台发电机中性点均不接地；二台变压器均为Y0/-11接线，发电机侧为接，Y0侧中性点直接接地；三条电力线路的零序电抗均
20、为0.20（以50MVA为基准值）。 要求计算节点3处a相接地短路，b、c两相短路，以及b、c两相接地短路时，故障点的短路电流。（Zf=Zg=0）五、源码1.用节点导纳矩阵计算三相短路电流clc;clear all;NF=1;n=4;nl=5;B=0 1 0.2i 0 1 0;0 2 4i 0 1 0; 1 3 0.51i 0 1 0; 2 3 0.59i 0 1 0; 3 4 1.43i 0 1 0;V0=1;1;1;1;D=4,0;m=0;Z=zeros(n);V=zeros(n);I=zeros(nl);for k1=1:nl p=B(k1,1); q=B(k1,2); if B(k1,
21、6)=0 %变压器支路首端是低压侧 l=1./B(k1,5); else l=B(k1,5);%变压器支路首端是高压侧 end if p=0 if qm %追加接地树支 Z(q,q)=B(k1,3); m=m+1; else %追加接地连支 for k=1:m, Z(k,m+1)=-Z(k,q); Z(m+1,k)=-Z(q,k); end Z(m+1,m+1)=Z(q,q)+B(k1,3); for l1=1:m for k=1:m Z(l1,k)=Z(l1,k)-Z(l1,m+1)*Z(m+1,k)./Z(m+1,m+1); end Z(l1,m+1)=0; end for k=1:m+1
22、 Z(m+1,k)=0; end end else if qm %追加不接地树支 for k=1:m Z(k,q)=Z(k,p)*l; Z(q,k)=Z(p,k)*l; end Z(q,q)=l2*Z(p,p)+l2*B(k1,3); m=m+1; else for k=1:m %追加不接地连支 Z(k,m+1)=l*Z(k,p)-Z(k,q); Z(m+1,k)=l*Z(p,k)-Z(q,k); end Z(m+1,m+1)=l2*Z(p,p)+Z(q,q)-2*l*Z(p,q)+l2*B(k1,3); for l1=1:m for k=1:m Z(l1,k)=Z(l1,k)-Z(l1,m+
23、1)*Z(m+1,k)./Z(m+1,m+1); end Z(l1,m+1)=0; end for k=1:m+1 Z(m+1,k)=0; end end endendfor k=1:NF I(D(k,1),D(k,1)=V0(D(k,1),1)./(Z(D(k,1),D(k,1)+D(k,2); ft=num2str(D(k,1); ts1=(点短路时); ts2=(电流的标幺值If=); dn=strcat(ft,ts1,ts2); disp(dn); disp( I(D(k,1),D(k,1); for i=1:n V(i,i)=V0(i,1)-I(D(k,1),D(k,1)*Z(i,D
24、(k,1); end for i=1:nl if B(i,6)=0 k=B(i,5); else k=1./B(i,5); end p=B(i,1);q=B(i,2); if p=0 e=0;b=B(i,3); I(i,i)=(e-V(q,q)./k)./b; else I(i,i)=(V(p,p)-V(q,q)./k)./B(i,3); end end disp(各节点的电压标幺值U为(节点号从小到大排):); for i=1:n disp(V(i,i); end disp(各支路短路电流的标幺值I为(顺序同您输入B时一样):) for i=1:nl disp(I(i,i); endend运
25、行结果：4点短路时电流的标幺值If= 0 - 0.4890i各节点的电压标幺值U为(节点号从小到大排): 0.9153 0.7380 0.6993 0各支路短路电流的标幺值I为(顺序同您输入B时一样): 0 + 4.5765i 0 + 0.1845i 0 - 0.4235i 0 - 0.0655i 0 - 0.4890i2.用计算曲线计算三相短路电流clc;clear all;NF=1;n=3;nl=5;Sb=100;D=3 0;vb=115 115 115;V0=1;1;1;B=0 1 0.15i 0 1 0 100 1 1;0 2 0.075i 0 1 0 200 1 1; 1 2 0.1
26、i 0 1 0 0 0 2; 1 3 0.1i 0 1 0 0 0 2; 2 3 0.1i 0 1 0 0 0 2;N=6;m=0;Z=zeros(n);sq3=sqrt(3);s=zeros(N,1);for k1=1:nl p=B(k1,1); q=B(k1,2); if B(k1,6)=0 %变压器支路首端是低压侧 l=1./B(k1,5); else l=B(k1,5);%变压器支路首端是高压侧 end if p=0 if qm %追加接地树支 Z(q,q)=B(k1,3); m=m+1; else %追加接地连支 for k=1:m, Z(k,m+1)=-Z(k,q); Z(m+1,
27、k)=-Z(q,k); end Z(m+1,m+1)=Z(q,q)+B(k1,3); for l1=1:m for k=1:m Z(l1,k)=Z(l1,k)-Z(l1,m+1)*Z(m+1,k)./Z(m+1,m+1); end Z(l1,m+1)=0; end for k=1:m+1 Z(m+1,k)=0; end end else if qm %追加不接地树支 for k=1:m Z(k,q)=Z(k,p)*l; Z(q,k)=Z(p,k)*l; end Z(q,q)=l2*Z(p,p)+l2*B(k1,3); m=m+1; else for k=1:m %追加不接地连支 Z(k,m+1)=l*Z(k,p)-Z(k,q); Z(m+1,k)=l*Z(p,k)-Z(q,k); end Z(m+1,m+1)=l2*Z(p,p)+Z(q,q)-2*l*Z(p,q)+l2*B(k1,3); for l1=1:m for k=1:m Z(l1,k)=Z(l1,k)-Z(l1,m+1)*Z(m+1,k)./Z(m+1,m+1); end Z(l1,m+1)=0; end for k=1:m+1 Z(m+1,k)=0; end end endendVb=zeros(1,n);V=zeros(1,n);Vd=zeros(1,n);
